Skip to content

Instantly share code, notes, and snippets.

@ntustison
Last active September 18, 2024 16:48
Show Gist options
  • Save ntustison/12a656a5fc2f6f9c4494c88dc09c5621 to your computer and use it in GitHub Desktop.
Save ntustison/12a656a5fc2f6f9c4494c88dc09c5621 to your computer and use it in GitHub Desktop.
ANTsX Tutorial

ANTs intro

  • ANTs --- Set of C++ programs built on top of the Insight Toolkit
  • ANTsR --- R-based ANTs interface (dependencies include ITKR, ANTsRCore)
  • ANTsPy --- Python-based ANTs interface
  • ANTsXNet --- deep learning library built on Keras
    • ANTsRNet
    • ANTsPyNet

ANTsR / ANTsPy

ANTsXNet


ANTsR / ANTsPy

Installation

  • Python: git clone https://github.com/ANTsX/ANTsPy; cd ANTsPy; python3 setup.py install
  • R: devtools::install_github( "ANTsX/ANTsR" )

Basics

Read, write, and visualize images

Python

>>> import ants
>>>
>>> r16_filename = ants.get_ants_data("r16")
>>> r16_image = ants.image_read(r16_filename)
>>>
>>> ants.write_image(r16_image, "r16.nii.gz")
>>>
>>> ants.plot( r16_image )

R

> library( ANTsR )
>
> r16Filename <- getANTsRData( "r16" )
> r16Image <- antsImageRead( r16Filename )
>
> antsImageWrite( r16Image, "r16.nii.gz" )
>
> plot( r16Image )

Image info

Python

>>> import ants
>>>
>>> r16 = ants.image_read(ants.get_ants_data('r16'))
>>> r16
ANTsImage
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (256, 256)
	 Spacing    : (1.0, 1.0)
	 Origin     : (0.0, 0.0)
	 Direction  : [1. 0. 0. 1.]
>>> r16.shape
(256, 256)
>>> r16.origin
(0.0, 0.0)
>>> ants.get_origin(r16)
(0.0, 0.0)
>>> r16.direction
array([[1., 0.],
       [0., 1.]])
>>> ants.get_direction(r16)
array([[1., 0.],
       [0., 1.]])
>>> r16.spacing
(1.0, 1.0)              
>>> ants.get_spacing(r16)
(1.0, 1.0)

R

> library( ANTsR )
>
> r16 <- antsImageRead( getANTsRData( 'r16' ) )
> r16
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 256x256 
  Voxel Spacing       : 1x1 
  Origin              : 0 0 
  Direction           : 1 0 0 1 
  Filename           : /Library/Frameworks/R.framework/Versions/4.0/Resources/library/ANTsRCore/extdata/r16slice.jpg 
> dim( r16 )
[1] 256 256
> antsGetOrigin( r16 )
[1] 0 0
> antsGetDirection( r16 )
     [,1] [,2]
[1,]    1    0
[2,]    0    1
> antsGetSpacing( r16 )
[1] 1 1

Conversion to/from native types

Python

>>> import ants
>>>
>>> r16 = ants.image_read(ants.get_ants_data('r16'))
>>> r16_numpy = r16.numpy()
>>> r16_numpy
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)
>>> r16_image = ants.from_numpy(r16_numpy, origin=r16.origin, spacing=r16.spacing, direction=r16.direction)       

R

> library( ANTsR )
>
> r16 <- antsImageRead( getANTsRData( 'r16' ) )
> r16Array <- as.array( r16 )
> r16Image <- as.antsImage( r16Array, spacing = antsGetSpacing( r16 ), origin = antsGetOrigin( r16 ), direction = antsGetDirection( r16 ) )
> r16Image <- as.antsImage( r16Array, reference = r16 )

Arithmetic operations

Python

>>> import ants
>>>
>>> r16 = ants.image_read(ants.get_ants_data('r16'))
>>> r64 = ants.image_read(ants.get_ants_data('r64'))
>>>
>>> r_x = r16 + r64
>>> r_x = r16 - r64
>>> r_x = r16 * r64
>>> r_x = r16 / r64
>>>
>>> r_x = r16 + 3.14
>>> r_x = r16 - 3.14
>>> r_x = r16 * 3.14
>>> r_x = r16 / 3.14
>>> r_x = r16 ** 3.14

R

> library( ANTsR )
>
> r16 <- antsImageRead( getANTsRData( "r16" ) )
> r64 <- antsImageRead( getANTsRData( "r64" ) )
>
> rx <- r16 + r64
> rx <- r16 - r64
> rx <- r16 * r64
> rx <- r16 / r64
>
> rx <- r16 + 3.14
> rx <- r16 - 3.14
> rx <- r16 * 3.14
> rx <- r16 / 3.14
> rx <- r16 ^ 3.14

Other operations

Python

>>> import ants
>>>
>>> r16 = ants.image_read(ants.get_ants_data('r16'))
>>>
>>> r16.min()
0.0
>>> r16.max()
254.0
>>> r16.std()
84.8173
>>> r16.mean()
50.461365

R

> library( ANTsR )
>
> r16 <- antsImageRead( getANTsRData( "r16" ) )
>
> min( r16 )
[1] 0
> max( r16 )
[1] 254
> sd( r16 )
[1] 84.81795
> mean( r16 )
[1] 50.46136

ANTsR / ANTsPy

Image registration

Python

>>> import ants
>>>
>>> r16 = ants.image_read(ants.get_ants_data('r16'))
>>> r16_seg = ants.threshold_image(r16, "Kmeans", 3)
 Kmeans with 3 thresholds
>>> ants.plot(r16, overlay=r16_seg, overlay_alpha=0.5)
>>>
>>> r64 = ants.image_read(ants.get_ants_data('r64'))
>>> r64_seg = ants.threshold_image(r64, "Kmeans", 3)
 Kmeans with 3 thresholds
>>> ants.plot(r64, overlay=r64_seg, overlay_alpha=0.5)
>>>
>>> # Affine only
>>> reg = ants.registration(fixed=r16, moving=r64, type_of_transform="antsRegistrationSyNQuick[a]", verbose=True)
>>> r16xr64_seg = ants.apply_transforms(fixed=r16, moving=r64_seg, transformlist=reg['fwdtransforms'], interpolator='nearestNeighbor')
>>> ants.plot(r16, overlay=r16xr64_seg, overlay_alpha=0.5)
>>> r64xr16_seg = ants.apply_transforms(fixed=r64, moving=r16_seg, transformlist=reg['invtransforms'], whichtoinvert=[True], interpolator='nearestNeighbor')
>>> ants.plot(r64, overlay=r64xr16_seg, overlay_alpha=0.5)
>>>
>>> # Linear and deformable
>>> reg = ants.registration(fixed=r16, moving=r64, type_of_transform="antsRegistrationSyNQuick[b]", verbose=True)
>>> r16xr64_seg = ants.apply_transforms(fixed=r16, moving=r64_seg, transformlist=reg['fwdtransforms'], interpolator='nearestNeighbor')
>>> ants.plot(r16, overlay=r16xr64_seg, overlay_alpha=0.5)
>>> r64xr16_seg = ants.apply_transforms(fixed=r64, moving=r16_seg, transformlist=reg['invtransforms'], whichtoinvert=[True, False], interpolator='nearestNeighbor')
>>> ants.plot(r64, overlay=r64xr16_seg, overlay_alpha=0.5)
>>>
>>> # Multivariate linear and deformable using intensity with wm and gm segmentations
>>> r16_gm = ants.threshold_image(r16_seg, 3, 3, 1, 0)
>>> r64_gm = ants.threshold_image(r64_seg, 3, 3, 1, 0)
>>> r16_wm = ants.threshold_image(r16_seg, 4, 4, 1, 0)
>>> r64_wm = ants.threshold_image(r64_seg, 4, 4, 1, 0)
>>> multivariate_extras = list()
>>> multivariate_extras.append(['MSQ', r16_gm, r64_gm, 2, 1])
>>> multivariate_extras.append(['MSQ', r16_wm, r64_wm, 2, 1])
>>> reg = ants.registration(fixed=r16, moving=r64, type_of_transform="antsRegistrationSyNQuick[b]", multivariate_extras = multivariate_extras, verbose=True)
>>> r16xr64_seg = ants.apply_transforms(fixed=r16, moving=r64_seg, transformlist=reg['fwdtransforms'], interpolator='nearestNeighbor')
>>> ants.plot(r16, overlay=r16xr64_seg, overlay_alpha=0.5)
>>> r64xr16_seg = ants.apply_transforms(fixed=r64, moving=r16_seg, transformlist=reg['invtransforms'], whichtoinvert=[True, False], interpolator='nearestNeighbor')
>>> ants.plot(r64, overlay=r64xr16_seg, overlay_alpha=0.5)

R

> library( ANTsR )
>
> r16 <- antsImageRead( getANTsRData( 'r16' ) )
> r16Seg <- thresholdImage( r16, "Kmeans", 3 )
 Kmeans with 3 thresholds
> plot( r16, r16Seg, alpha = 0.5 )
>
> r64 <- antsImageRead( getANTsRData( 'r64' ) )
> r64Seg <- thresholdImage( r64, "Kmeans", 3 )
 Kmeans with 3 thresholds
> plot( r64, r64Seg, alpha = 0.5 )
>
> # Affine only
> reg <- antsRegistration( fixed = r16, moving = r64, typeofTransform = "antsRegistrationSyNQuick[a]", verbose = TRUE )
> r16xr64Seg <- antsApplyTransforms( fixed = r16, moving = r64Seg, transformlist = reg$fwdtransforms, interpolator = 'nearestNeighbor' )
> plot( r16, r16xr64Seg, alpha = 0.5 ) 
> r64xr16Seg <- antsApplyTransforms( fixed = r64, moving = r16Seg, transformlist = reg$invtransforms, whichtoinvert = c( TRUE ), interpolator = 'nearestNeighbor' )
> plot( r64, r64xr16Seg, alpha = 0.5 ) 
>
> # Linear and deformable
> reg <- antsRegistration( fixed = r16, moving = r64, typeofTransform = "antsRegistrationSyNQuick[b]", verbose = TRUE )
> r16xr64Seg <- antsApplyTransforms( fixed = r16, moving = r64Seg, transformlist = reg$fwdtransforms, interpolator = 'nearestNeighbor' )
> plot( r16, r16xr64Seg, alpha = 0.5 ) 
> r64xr16Seg <- antsApplyTransforms( fixed = r64, moving = r16Seg, transformlist = reg$invtransforms, whichtoinvert = c( TRUE, FALSE ), interpolator = 'nearestNeighbor' )
> plot( r64, r64xr16Seg, alpha = 0.5 ) 
>
> # Linear and deformable using intensity with wm and gm segmentations
> r16Gm <- thresholdImage( r16Seg, 3, 3, 1, 0 )
> r64Gm <- thresholdImage( r64Seg, 3, 3, 1, 0 )
> r16Wm <- thresholdImage( r16Seg, 4, 4, 1, 0 )
> r64Wm <- thresholdImage( r64Seg, 4, 4, 1, 0 ) 
> multivariateExtras <- list()
> multivariateExtras[[1]] <- list( 'MSQ', r16Gm, r64Gm, 2, 1 )
> multivariateExtras[[2]] <- list( 'MSQ', r16Wm, r64Wm, 2, 1 )
> reg <- antsRegistration( fixed = r16, moving = r64, typeofTransform = "antsRegistrationSyNQuick[b]", multivariateExtras = multivariateExtras, verbose = TRUE )
> r16xr64Seg <- antsApplyTransforms( fixed = r16, moving = r64Seg, transformlist = reg$fwdtransforms, interpolator = 'nearestNeighbor' )
> plot( r16, r16xr64Seg, alpha = 0.5 ) 
> r64xr16Seg <- antsApplyTransforms( fixed = r64, moving = r16Seg, transformlist = reg$invtransforms, whichtoinvert = c( TRUE, FALSE ), interpolator = 'nearestNeighbor' )
> plot( r64, r64xr16Seg, alpha = 0.5 ) 

Template building

Python

>>> import ants
>>>
>>> r16 = ants.image_read(ants.get_ants_data('r16'))
>>> r27 = ants.image_read(ants.get_ants_data('r27'))
>>> r30 = ants.image_read(ants.get_ants_data('r30'))
>>> r62 = ants.image_read(ants.get_ants_data('r62'))
>>> r64 = ants.image_read(ants.get_ants_data('r64'))
>>> r85 = ants.image_read(ants.get_ants_data('r85'))
>>>
>>> rimage_list = [r16, r27, r30, r62, r64, r85]
>>> rtemplate = ants.build_template(initial_template=None, image_list=rimage_list, type_of_transform="antsRegistrationSyNQuick[s]")
>>> ants.plot(rtemplate)

R

> library( ANTsR )
> 
> r16 <- antsImageRead( getANTsRData( 'r16' ) )
> r27 <- antsImageRead( getANTsRData( 'r27' ) )
> r30 <- antsImageRead( getANTsRData( 'r30' ) )
> r62 <- antsImageRead( getANTsRData( 'r62' ) )
> r64 <- antsImageRead( getANTsRData( 'r64' ) )
> r85 <- antsImageRead( getANTsRData( 'r85' ) )
> 
> rimageList <- list( r16, r27, r30, r62, r64, r85 )
> initialAverage <- antsAverageImages( rimageList )
> rtemplate <- buildTemplate( initialTemplate = initialAverage, imgList = rimageList, typeofTransform = "antsRegistrationSyNQuick[s]", verbose = TRUE )
> plot( rtemplate )

Motion correction

Python

>>> import ants
>>> 
>>> # Reduce the number of times points
>>> image = ants.image_read(ants.get_ants_data("pcasl"))
>>> image_array = image.numpy()
>>> image_array = image_array[:,:,:,:10]
>>> image = ants.from_numpy(image_array, spacing=image.spacing, direction=image.direction, origin=image.origin)
>>> 
>>> moco = ants.motion_correction(image, type_of_transform="BOLDRigid", verbose=True)
>>> moco['motion_corrected']
ANTsImage
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (64, 64, 18, 10)
	 Spacing    : (3.4375, 3.4375, 7.2, 4.0)
	 Origin     : (-116.0013, 57.1198, -38.3576, 0.0)
	 Direction  : [ 1.  0.  0.  0.  0. -1.  0.  0.  0.  0.  1.  0.  0.  0.  0.  1.]
>>> moco['motion_parameters']
[['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmpy49zfa890GenericAffine.mat'], ['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmp_qncwnum0GenericAffine.mat'], ['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmpwlnjbm_30GenericAffine.mat'], ['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmppfj9lmyv0GenericAffine.mat'], ['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmpw5y0z0ot0GenericAffine.mat'], ['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmpztsq1ete0GenericAffine.mat'], ['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmptw81m6lx0GenericAffine.mat'], ['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmptxnjmkc_0GenericAffine.mat'], ['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmpwuv83l_a0GenericAffine.mat'], ['/var/folders/5h/62f7bql14zvgx34nqwd7xvnc0000gn/T/tmps8_zjfy60GenericAffine.mat']]
>>> moco['FD']
array([0.        , 0.09816545, 0.18367779, 0.12037325, 0.13567556,
       0.09293036, 0.09830562, 0.18831979, 0.04102131, 0.06663223])

R

> library( ANTsR )
> 
> # Reduce the number of times points
> image <- antsImageRead( getANTsRData( "pcasl" ) )
> imageArray <- as.array( image )
> imageArray <- imageArray[,,,1:10]
> image <- as.antsImage( imageArray, reference = image )
> 
> moco <- antsMotionCalculation( image, txtype = "Rigid", verbose = TRUE )
> moco$moco_img
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 64x64x18x10 
  Voxel Spacing       : 3.4375x3.4375x7.19999980926514x4 
  Origin              : -116.0013 57.11981 -38.35765 0 
  Direction           : 1 0 0 0 0 -1 0 0 0 0 1 0 0 0 0 1 
> moco$moco_params
   MetricPre MetricPost     MOCOparam1     MOCOparam2     MOCOparam3
1          0  -1.005069  0.00023146855 -0.00009585432  0.00004607494
2          0  -1.005325 -0.00107804559 -0.00067556275  0.00025181432
3          0  -1.009195 -0.00018273871  0.00001386437 -0.00002405675
4          0  -1.010406 -0.00015799195 -0.00019400520  0.00046687210
5          0  -1.007293  0.00011430427 -0.00004679069 -0.00001049459
6          0  -1.005072  0.00022220720 -0.00007070929 -0.00003112626
7          0  -1.008751  0.00009608059 -0.00001348405  0.00001652053
8          0  -1.003125 -0.00132131224  0.00013327676 -0.00000959006
9          0  -1.005963 -0.00009936649  0.00001233287 -0.00002474547
10         0  -1.007560 -0.00061640351 -0.00064586266 -0.00052936993
     MOCOparam4   MOCOparam5   MOCOparam6
1   0.010914629 -0.008945767 -0.111323305
2   0.001277686 -0.079050128 -0.091623102
3  -0.021036213  0.014788211  0.098308572
4  -0.003656350  0.037733644  0.004530848
5   0.007828627  0.008543149 -0.075201220
6  -0.018174798  0.003812138 -0.146378820
7   0.015950238 -0.030493634 -0.063736015
8  -0.010916564  0.004372850 -0.006914493
9  -0.002273881  0.022307026  0.068820824
10  0.027577054 -0.036157711  0.055110048
> moco$fd
   MeanDisplacement MaxDisplacement
1        0.09661564       0.1930724
2        0.18165174       0.2670051
3        0.10522685       0.1257340
4        0.09715838       0.1322308
5        0.08049656       0.0901812
6        0.10052631       0.1125070
7        0.13801552       0.2513173
8        0.05978751       0.1279247
9        0.06027606       0.1109254
10       0.00000000       0.0000000

Landmark registration

Python

Two simple examples

>>> import ants
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> 
>>> examples = ("Big square to little square",
>>>             "Brian's canonical example")
>>> 
>>> moving_points_list = list()
>>> fixed_points_list = list()
>>> 
>>> # Big square to little square
>>> moving_points = np.zeros((4, 2))
>>> moving_points[0, 0] = 64
>>> moving_points[0, 1] = 64
>>> moving_points[1, 0] = 192
>>> moving_points[1, 1] = 64
>>> moving_points[2, 0] = 64
>>> moving_points[2, 1] = 192
>>> moving_points[3, 0] = 192
>>> moving_points[3, 1] = 192
>>> moving_points_list.append(moving_points)
>>> 
>>> fixed_points = np.zeros((4, 2))
>>> fixed_points[0, 0] = 96
>>> fixed_points[0, 1] = 96
>>> fixed_points[1, 0] = 160
>>> fixed_points[1, 1] = 96
>>> fixed_points[2, 0] = 96
>>> fixed_points[2, 1] = 160
>>> fixed_points[3, 0] = 160
>>> fixed_points[3, 1] = 160
>>> fixed_points_list.append(fixed_points)
>>> 
>>> # Brian's example
>>> moving_points = np.zeros((3, 2))
>>> moving_points[0, 0] = 64
>>> moving_points[0, 1] = 64
>>> moving_points[1, 0] = 192
>>> moving_points[1, 1] = 64
>>> moving_points[2, 0] = 64
>>> moving_points[2, 1] = 192
>>> moving_points_list.append(moving_points)
>>> 
>>> fixed_points = np.zeros((3, 2))
>>> fixed_points[0, 0] = 192
>>> fixed_points[0, 1] = 192
>>> fixed_points[1, 0] = 192
>>> fixed_points[1, 1] = 64
>>> fixed_points[2, 0] = 64
>>> fixed_points[2, 1] = 192
>>> fixed_points_list.append(fixed_points)
>>> 
>>> # domain image
>>> r16 = ants.image_read(ants.get_ants_data("r16"))
>>> 
>>> transform_types = ('rigid', 'similarity', 'affine', 'bspline', 'tps', 'diffeo', 'syn', 'time-varying')
>>> 
>>> for e in range(len(examples)):
>>>     warped_grid_arrays = list()
>>>     for i in range(len(transform_types)):
>>>         print("Fitting with transform ", transform_types[i])
>>>         xfrm = ants.fit_transform_to_paired_points(moving_points_list[e], fixed_points_list[e], 
>>>                            transform_type=transform_types[i],
>>>                           domain_image=r16, 
>>>                           number_of_fitting_levels=4, 
>>>                           mesh_size=(4, 4), 
>>>                           number_of_compositions=10000,
>>>                           number_of_integration_steps=10,
>>>                           number_of_time_steps=2,
>>>                           composition_step_size=0.2,
>>>                           convergence_threshold=1e-6,
>>>                           enforce_stationary_boundary=True, 
>>>                           verbose=True)
>>>        grid = ants.create_warped_grid(r16)
>>>        if isinstance(xfrm, dict):
>>>            warped_grid = xfrm['forward_transform'].apply_to_image(grid, r16)
>>>        else:    
>>>            warped_grid = xfrm.apply_to_image(grid, r16)
>>>        warped_grid_arrays.append(warped_grid.numpy()) 
>>>    fig, axes = plt.subplots(2, 4, figsize=(10, 10))
>>>    for i, ax in enumerate(axes.flat):
>>>        ax.imshow(warped_grid_arrays[i], cmap="gray")
>>>        if transform_types[i] == "bspline":
>>>            ax.set_title("Transform: bspline (clamped)")
>>>        else:
>>>            ax.set_title("Transform: " + transform_types[i])
>>>        ax.xaxis.set_visible(False)
>>>        ax.yaxis.set_visible(False)
>>>        for j in range(moving_points_list[e].shape[0]):
>>>            ax.arrow(moving_points_list[e][j,0], moving_points_list[e][j,1],
>>>                     fixed_points_list[e][j,0]-moving_points_list[e][j,0],
>>>                     fixed_points_list[e][j,1]-moving_points_list[e][j,1],
>>>                     head_width=10.0,
>>>                     width=2.0,
>>>                     color='red',
>>>                     ec='red')
>>>
>>>    fig.suptitle(examples[e])
>>>    fig.tight_layout()
>>>    plt.show()
>>>    plt.close()

More complicated example

>>> import ants
>>> import numpy as np
>>> import os
>>> import matplotlib.pyplot as plt
>>> 
>>> os.environ["ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS"] = "2"
>>> 
>>> x = np.linspace(-1, 1, num=20)
>>> 
>>> moving_x = np.empty(x.shape)
>>> moving_y = np.empty(x.shape)
>>> fixed_x = np.empty(x.shape)
>>> fixed_y = np.empty(x.shape)
>>> 
>>> center_x = 128
>>> center_y = 128
>>> offset_x = 50
>>> offset_y = 50
>>> scale = 64
>>> 
>>> for i in range(len(x)):
>>>     fixed_x[i] = x[i] * scale + center_x + offset_x
>>>     fixed_y[i] = x[i] ** 2 * scale + center_y + offset_y
>>>     moving_x[i] = x[i] * scale + center_x
>>>     moving_y[i] = -x[i] ** 2 * scale + center_y
>>> 
>>> fixed_points = np.vstack((fixed_x, fixed_y)).transpose()
>>> moving_points = np.vstack((moving_x, moving_y)).transpose()
>>> 
>>> #########   rigid transform    ############
>>> 
>>> rigid_xfrm = ants.fit_transform_to_paired_points(moving_points, fixed_points, transform_type="rigid")
>>> 
>>> fixed_points_rigid_warped = np.empty(fixed_points.shape)
>>> for j in range(fixed_points.shape[0]):
>>>     fixed_points_rigid_warped[j,:] = rigid_xfrm.apply_to_point(tuple(fixed_points[j,:]))
>>> 
>>> #########   syn transform    ############
>>> 
>>> r16 = ants.image_read(ants.get_ants_data("r16"))
>>> 
>>> syn_xfrm = ants.fit_transform_to_paired_points(moving_points, fixed_points_rigid_warped, 
>>>                       transform_type="syn",
>>>                       domain_image=r16, 
>>>                       number_of_fitting_levels=4, 
>>>                       mesh_size=(9, 9), 
>>>                       number_of_compositions=100,
>>>                       composition_step_size=0.2,
>>>                       enforce_stationary_boundary=True, 
>>>                       verbose=True)
>>> 
>>> total_forward_xfrm = ants.compose_ants_transforms([rigid_xfrm, syn_xfrm['forward_transform']])
>>> total_fixed_to_middle_xfrm = ants.compose_ants_transforms([rigid_xfrm, syn_xfrm['fixed_to_middle_transform']])
>>> 
>>> fixed_points_syn_warped = np.empty(fixed_points.shape)
>>> fixed_points_syn_mid_warped = np.empty(fixed_points.shape)
>>> moving_points_syn_warped = np.empty(fixed_points.shape)
>>> moving_points_syn_mid_warped = np.empty(fixed_points.shape)
>>> 
>>> for j in range(fixed_points.shape[0]):
>>>     fixed_points_syn_warped[j,:] = total_forward_xfrm.apply_to_point(tuple(fixed_points[j,:]))
>>>     fixed_points_syn_mid_warped[j,:] = total_fixed_to_middle_xfrm.apply_to_point(tuple(fixed_points[j,:]))
>>>     moving_points_syn_warped[j,:] = syn_xfrm['inverse_transform'].apply_to_point(tuple(moving_points[j,:]))
>>>     moving_points_syn_mid_warped[j,:] = syn_xfrm['moving_to_middle_transform'].apply_to_point(tuple(moving_points[j,:]))
>>> 
>>> plt.plot(fixed_points[:,0], fixed_points[:,1], 'g^', label="Fixed points") 
>>> plt.plot(moving_points[:,0], moving_points[:,1], 'go', label="Moving points")
>>> plt.plot(fixed_points_rigid_warped[:,0], fixed_points_rigid_warped[:,1], 'r^', label="Fixed (rigid)")
>>> plt.plot(fixed_points_syn_warped[:,0], fixed_points_syn_warped[:,1], 'r^--', alpha=0.33, label="Fixed (SyN-forward)") 
>>> plt.plot(moving_points_syn_warped[:,0], moving_points_syn_warped[:,1], 'ro--', alpha=0.33, label="Moving (SyN-inverse)") 
>>> plt.plot(fixed_points_syn_mid_warped[:,0], fixed_points_syn_mid_warped[:,1], 'b^', alpha=1.0, label="Fixed (SyN-mid)") 
>>> plt.plot(moving_points_syn_mid_warped[:,0], moving_points_syn_mid_warped[:,1], 'bo--', alpha=0.33, label="Moving (SyN-mid)")
>>> 
>>> plt.legend()
>>> plt.title('Rigid + SyN landmark registration')
>>> plt.show()
>>>
>>> grid = ants.create_warped_grid(r16)
>>> warped_grid = syn_xfrm['forward_transform'].apply_to_image(grid, r16)
>>> ants.plot(warped_grid)

R

Two simple examples

> library( ANTsR )
> library( ggplot2 )
> library( raster )
> 
> examples <- c("Big square to little square",
>               "Brian's canonical example")
> 
> movingPointsList <- list()
> fixedPointsList <- list()
> 
> # Big square to little square
> movingPoints <- array( data = 0, dim = c( 4, 2 ) )
> movingPoints[1, 1] <- 64
> movingPoints[1, 2] <- 64
> movingPoints[2, 1] <- 192
> movingPoints[2, 2] <- 64
> movingPoints[3, 1] <- 64
> movingPoints[3, 2] <- 192
> movingPoints[4, 1] <- 192
> movingPoints[4, 2] <- 192
> movingPointsList[[1]] <- movingPoints
> 
> fixedPoints <- array( data = 0, dim = c( 4, 2 ) )
> fixedPoints[1, 1] <- 96
> fixedPoints[1, 2] <- 96
> fixedPoints[2, 1] <- 160
> fixedPoints[2, 2] <- 96
> fixedPoints[3, 1] <- 96
> fixedPoints[3, 2] <- 160
> fixedPoints[4, 1] <- 160
> fixedPoints[4, 2] <- 160
> fixedPointsList[[1]] <- fixedPoints
> 
> # Brian's example
> movingPoints <- array( data = 0, dim = c( 3, 2 ) )
> movingPoints[1, 1] <- 64
> movingPoints[2, 1] <- 192
> movingPoints[2, 2] <- 64
> movingPoints[3, 1] <- 64
> movingPoints[3, 2] <- 192
> movingPointsList[[2]] <- movingPoints
>
> fixedPoints[1, 1] <- 192
> fixedPoints[1, 2] <- 192
> fixedPoints[2, 1] <- 192
> fixedPoints[2, 2] <- 64
> fixedPoints[3, 1] <- 64
> fixedPoints[3, 2] <- 192
> fixedPointsList[[2]] <- fixedPoints
> 
> # domain image
> r16 <- antsImageRead( getANTsRData( "r16" ) )
> 
> transformTypes <- c( 'rigid', 'similarity', 'affine', 'bspline', 'tps', 'diffeo', 'syn', 'time-varying' )
> 
> for( e in seq.int( length( examples ) ) )
>   {
>   warpedGridArrays <- list()
>   for( i in seq.int( length( transformTypes ) ) )
>     {
>     cat( paste0( "Fitting with transform ", transformTypes[i], "\n" ) )
>     xfrm <- fitTransformToPairedPoints( movingPointsList[[e]], fixedPointsList[[e]], 
>                        transformType = transformTypes[i],
>                        domainImage = r16, 
>                        numberOfFittingLevels = 4, 
>                        meshSize = c( 4, 4 ), 
>                        numberOfCompositions = 10000,
>                        numberOfIntegrationSteps = 10,
>                        numberOfTimeSteps = 2,
>                        compositionStepSize = 0.2,
>                        convergenceThreshold = 1e-6,
>                        enforceStationaryBoundary = TRUE, 
>                        verbose = TRUE )
>     grid <- createWarpedGrid( r16 )
>     if( is.list( xfrm ) )
>       {
>       warpedGrid <- applyAntsrTransformToImage( xfrm$forwardTransform, grid, r16 )
>       } else {
>       warpedGrid <- applyAntsrTransformToImage( xfrm, grid, r16 )
>       }
>     warpedGridArrays[[i]] <- as.array( warpedGrid )
>     }
> 
>   plotDataFrame <- data.frame()
>   for( i in seq.int( length( warpedGridArrays ) ) )  
>     {
>     rasterDataFrame <- as.data.frame( raster( warpedGridArrays[[i]] ), xy = TRUE )
>     rasterDataFrame$TransformType <- rep( transformTypes[i], nrow( rasterDataFrame ) )
>     if( i == 1 )
>       {
>       plotDataFrame <- rasterDataFrame
>       } else {
>       plotDataFrame <- rbind( plotDataFrame, rasterDataFrame )
>       }
>     }
> 
>   plotDataFrame$TransformType[which( plotDataFrame$TransformType == "bspline" )] <- "bspline (clamped)"
>   plotDataFrame$TransformType <- factor( plotDataFrame$TransformType, levels = unique( plotDataFrame$TransformType ) )  
> 
>   p <- ggplot( data = plotDataFrame, aes( x = x, y = y, fill = layer ), colour = "gray" ) + 
>        geom_raster() + 
>        facet_wrap( ~TransformType ) + 
>        scale_fill_gradientn( colours = c( "black", "white" ) ) +
>        theme( legend.position = "none",
>               axis.title.x = element_blank(), 
>               axis.title.y = element_blank(), 
>               axis.text.x = element_blank(), 
>               axis.text.y = element_blank(), 
>               axis.ticks.x = element_blank(), 
>               axis.ticks.y = element_blank()
>               ) +
>        ggtitle( examples[e] )      
> 
>   # ggsave( paste0( "~/Desktop/foo", e, ".pdf" ), plot = p )            
>   }

More complicated example

> library( ANTsR )
> library( ggplot2 )
> 
> Sys.setenv( ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS = "2" )
> 
> x <- seq( -1, 1, length.out = 20 )
> 
> movingX <- c()
> movingY <- c()
> fixedX <- c()
> fixedY <- c()
> 
> centerX <- 128
> centerY <- 128
> offsetY <- 50
> offsetX <- 50
> offsetY <- 50
> scale <- 64
> 
> for( i in seq.int( length( x ) ) )
>   {
>   fixedX[i] <- x[i] * scale + centerX + offsetX 
>   fixedY[i] <- x[i]^2 * scale + centerY + offsetY
>   movingX[i] <- x[i] * scale + centerX
>   movingY[i] <- -x[i]^2 * scale + centerY
>   }
> 
> fixedPoints <- cbind( fixedX, fixedY )
> movingPoints <- cbind( movingX, movingY )
> 
> #########   rigid transform    ############
> 
> rigidXfrm <- fitTransformToPairedPoints( movingPoints, fixedPoints, transformType = "rigid" )
> 
> fixedPointsRigidWarped <- applyAntsrTransform( rigidXfrm, fixedPoints )
> 
> #########   syn transform    ############
> 
> r16 <- antsImageRead( getANTsRData( "r16" ) )
> synXfrm <- fitTransformToPairedPoints( movingPoints, fixedPointsRigidWarped,
>                 transformType = "syn",
>                 domainImage = r16,
>                 numberOfFittingLevels = 4,
>                 meshSize = c( 9, 9 ),
>                 numberOfCompositions = 100,
>                 compositionStepSize = 0.2,
>                 enforceStationaryBoundary = TRUE,
>                 verbose = TRUE )
> 
> 
> totalForwardXfrm <- composeAntsrTransforms( list( rigidXfrm, synXfrm$forwardTransform ) )
> totalFixedToMiddleXfrm <- composeAntsrTransforms( list( rigidXfrm, synXfrm$fixedToMiddleTransform ) )
> 
> fixedPointsSynWarped <- applyAntsrTransform( totalForwardXfrm, fixedPoints )
> fixedPointsSynMidWarped <- applyAntsrTransform( totalFixedToMiddleXfrm, fixedPoints )
> movingPointsSynWarped <- applyAntsrTransform( synXfrm$inverseTransform, movingPoints )
> movingPointsSynMidWarped <- applyAntsrTransform( synXfrm$movingToMiddleTransform, movingPoints )
> 
> numberOfPoints <- dim( fixedPoints )[1]
> allPointsX <- c( fixedPoints[,1], movingPoints[,1], fixedPointsRigidWarped[,1],
>                  fixedPointsSynWarped[,1], movingPointsSynWarped[,1],
>                  fixedPointsSynMidWarped[,1], movingPointsSynMidWarped[,1] )
> allPointsY <- c( fixedPoints[,2], movingPoints[,2], fixedPointsRigidWarped[,2],
>                  fixedPointsSynWarped[,2], movingPointsSynWarped[,2],
>                  fixedPointsSynMidWarped[,2], movingPointsSynMidWarped[,2] )
> pointsType <- c( rep( "Fixed points", numberOfPoints ), rep( "Moving points", numberOfPoints ), rep( "Fixed (rigid)", numberOfPoints ), 
>                  rep( "Fixed (SyN-forward)", numberOfPoints ), rep( "Moving (SyN-inverse))", numberOfPoints ),
>                  rep( "Fixed (SyN-mid)", numberOfPoints ), rep( "Moving (SyN-mid))", numberOfPoints )
>                )
> 
> pointPlotDf <- data.frame( X = allPointsX, Y = allPointsY, Type = pointsType )
> pointPlot <- ggplot( data = pointPlotDf ) + 
>              geom_point( aes( x = X, y = Y, colour = Type, shape = Type ) ) + 
>              geom_line( aes( x = X, y = Y, colour = Type, group = Type ), alpha = 0.25 )
> pointPlot
> 
> grid <- createWarpedGrid( r16 )
> warpedGrid <- applyAntsrTransform( synXfrm$forwardTransform, grid )
> plot( warpedGrid )

Velocity flows across point sets

A more complicated example.

Python

>>> import ants
>>> 
>>> import math
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> 
>>> domain_array = np.zeros((100, 100))
>>> domain_image = ants.from_numpy(domain_array)
>>> 
>>> number_of_points = 8
>>> dimensionality = 2
>>> 
>>> rectangle = np.zeros((number_of_points, dimensionality))
>>> rectangle[0, :] = [35.0, 15.0]
>>> rectangle[1, :] = [35.0, 50.0]
>>> rectangle[2, :] = [35.0, 85.0]
>>> rectangle[3, :] = [50.0, 85.0]
>>> rectangle[4, :] = [65.0, 85.0]
>>> rectangle[5, :] = [65.0, 50.0]
>>> rectangle[6, :] = [65.0, 15.0]
>>> rectangle[7, :] = [50.0, 15.0]
>>> 
>>> square = np.zeros((number_of_points, dimensionality))
>>> square[0, :] = [25.0, 25.0]
>>> square[1, :] = [25.0, 50.0]
>>> square[2, :] = [25.0, 75.0]
>>> square[3, :] = [50.0, 75.0]
>>> square[4, :] = [75.0, 75.0]
>>> square[5, :] = [75.0, 50.0]
>>> square[6, :] = [75.0, 25.0]
>>> square[7, :] = [50.0, 25.0]
>>> 
>>> circle = np.zeros((number_of_points, dimensionality))
>>> for i in range(number_of_points):
>>>     circle[i, 0] = 50.0 + 25.0 * math.sin(2 * math.pi * i / 8 - 3 * math.pi / 4)
>>>     circle[i, 1] = 50.0 + 25.0 * math.cos(2 * math.pi * i / 8 - 3 * math.pi / 4)
>>> 
>>> point_sets = list()
>>> point_sets.append(rectangle)
>>> point_sets.append(square)
>>> point_sets.append(circle)
>>> 
>>> plt.plot(rectangle[:,0], rectangle[:,1], 'r-', label="Rectangle")
>>> plt.plot(square[:,0], square[:,1], 'g-', label="Square")
>>> plt.plot(circle[:,0], circle[:,1], 'b-', label="Circle")
>>> 
>>> plt.legend()
>>> plt.title('Original data')
>>> plt.show()
>>> 
>>> tv = ants.fit_time_varying_transform_to_point_sets(point_sets, time_points=None,
>>>   number_of_time_steps=None, domain_image=domain_image,
>>>   number_of_fitting_levels=4, mesh_size=2, spline_order=3, displacement_weights=None,
>>>   number_of_compositions=100, composition_step_size=0.5, number_of_integration_steps=10,
>>>   sigma=0.0, convergence_threshold=0.0, rasterize_points=False, verbose=True)
>>> 
>>> # Warp between the endpoints
>>> 
>>> warped_rectangle_to_circle = np.empty_like(rectangle)
>>> warped_rectangle_to_circle[:] = rectangle
>>> for j in range(rectangle.shape[0]):
>>>     warped_rectangle_to_circle[j,:] = tv['forward_transform'].apply_to_point(tuple(rectangle[j,:]))
>>> 
>>> warped_circle_to_rectangle = np.empty_like(circle)
>>> warped_circle_to_rectangle[:] = circle
>>> for j in range(circle.shape[0]):
>>>     warped_circle_to_rectangle[j,:] = tv['inverse_transform'].apply_to_point(tuple(circle[j,:]))
>>> 
>>> plt.plot(warped_rectangle_to_circle[:,0], warped_rectangle_to_circle[:,1], 'r--', label="RectangleToCircle")
>>> plt.plot(circle[:,0], circle[:,1], 'r-', label="Circle")
>>> plt.plot(warped_circle_to_rectangle[:,0], warped_circle_to_rectangle[:,1], 'b--', label="CircleToRectangle")
>>> plt.plot(rectangle[:,0], rectangle[:,1], 'b-', label="Rectangle")
>>> 
>>> plt.legend()
>>> plt.title('Warping between endpoints')
>>> plt.show()
>>> 
>>> # Warp the endpoints (i.e., rectangle and circle) to the middle (i.e., square)
>>> 
>>> forward_xfrm = ants.transform_from_displacement_field(ants.integrate_velocity_field(tv['velocity_field'], 0, 0.5))
>>> warped_rectangle_to_square = np.empty_like(rectangle)
>>> warped_rectangle_to_square[:] = rectangle
>>> for j in range(rectangle.shape[0]):
>>>     warped_rectangle_to_square[j,:] = forward_xfrm.apply_to_point(tuple(rectangle[j,:]))
>>> 
>>> inverse_xfrm = ants.transform_from_displacement_field(ants.integrate_velocity_field(tv['velocity_field'], 1.0, 0.5))
>>> warped_circle_to_square = np.empty_like(circle)
>>> warped_circle_to_square[:] = circle
>>> for j in range(circle.shape[0]):
>>>     warped_circle_to_square[j,:] = inverse_xfrm.apply_to_point(tuple(circle[j,:]))
>>> 
>>> plt.plot(warped_rectangle_to_square[:,0], warped_rectangle_to_square[:,1], 'r--', label="RectangleToSquare")
>>> plt.plot(warped_circle_to_square[:,0], warped_circle_to_square[:,1], 'b--', label="CircleToSquare")
>>> plt.plot(rectangle[:,0], rectangle[:,1], 'r-', label="Rectangle")
>>> plt.plot(square[:,0], square[:,1], 'g-', label="Square")
>>> plt.plot(circle[:,0], circle[:,1], 'b-', label="Circle")
>>> 
>>> plt.legend()
>>> plt.title('Warping to middle')
>>> plt.show()

R

> library( ANTsR )
> library( ggplot2 )
> 
> domainArray <- array( data = 0, dim = c( 100, 100 ) )
> domainImage <- as.antsImage( domainArray )
> 
> numberOfPoints <- 8
> dimensionality <- 2
> 
> rectangle <- array( data = 0, dim = c( numberOfPoints, dimensionality ) )
> rectangle[1,] <- c( 35.0, 15.0 )
> rectangle[2,] <- c( 35.0, 50.0 )
> rectangle[3,] <- c( 35.0, 85.0 )
> rectangle[4,] <- c( 50.0, 85.0 )
> rectangle[5,] <- c( 65.0, 85.0 )
> rectangle[6,] <- c( 65.0, 50.0 )
> rectangle[7,] <- c( 65.0, 15.0 )
> rectangle[8,] <- c( 50.0, 15.0 )
> 
> square <- array( data = 0, dim = c( numberOfPoints, dimensionality ) )
> square[1,] <- c( 25.0, 25.0 )
> square[2,] <- c( 25.0, 50.0 )
> square[3,] <- c( 25.0, 75.0 )
> square[4,] <- c( 50.0, 75.0 )
> square[5,] <- c( 75.0, 75.0 )
> square[6,] <- c( 75.0, 50.0 )
> square[7,] <- c( 75.0, 25.0 )
> square[8,] <- c( 50.0, 25.0 )
> 
> circle <- array( data = 0, dim = c( numberOfPoints, dimensionality ) )
> for( i in seq_len( numberOfPoints ) )
>   {
>   circle[i, 1] <- 50.0 + 25.0 * sin( 2 * pi * ( i - 1 ) / 8 - 3 * pi / 4 )
>   circle[i, 2] <- 50.0 + 25.0 * cos( 2 * pi * ( i - 1 ) / 8 - 3 * pi / 4 )
>   }
> 
> pointSets <- list()
> pointSets[[1]] <- rectangle
> pointSets[[2]] <- square
> pointSets[[3]] <- circle
> 
> allPointsX <- c( rectangle[,1], square[,1], circle[,1] )
> allPointsY <- c( rectangle[,2], square[,2], circle[,2] )
> pointsType <- c( rep( "Rectangle", numberOfPoints ), rep( "Square", numberOfPoints ), rep( "Circle", numberOfPoints ) )
> 
> pointPlotDf <- data.frame( X = allPointsX, Y = allPointsY, Type = pointsType )
> pointPlot <- ggplot( data = pointPlotDf ) +
>              geom_point( aes( x = X, y = Y, colour = Type, shape = Type ) ) +
>              geom_path( aes( x = X, y = Y, colour = Type ), alpha = 0.5 ) +
>              ggtitle( "Original data" )
> 
> pointPlot
> 
> tv <- fitTimeVaryingTransformToPointSets( pointSets, timePoints = NULL,
>   numberOfTimeSteps = 3, domainImage = domainImage,
>   numberOfFittingLevels = 4, meshSize = 2, splineOrder = 3, displacementWeights = NULL,
>   numberOfCompositions = 100, compositionStepSize = 0.5, numberOfIntegrationSteps = 10,
>   sigma = 0.0, convergenceThreshold = 0.0, rasterizePoints = FALSE, verbose = TRUE )
> 
> # Warp between the endpoints
> 
> warpedRectangleToCircle <- applyAntsrTransformToPoint( tv$forwardTransform, rectangle )
> warpedCircleToRectangle <- applyAntsrTransformToPoint( tv$inverseTransform, circle )
> 
> allPointsX <- c( warpedRectangleToCircle[,1], circle[,1], warpedCircleToRectangle[,1], rectangle[,1] )
> allPointsY <- c( warpedRectangleToCircle[,2], circle[,2], warpedCircleToRectangle[,2], rectangle[,2] )
> pointsType <- c( rep( "RectangleToCircle", numberOfPoints ), rep( "Circle", numberOfPoints ),
>                  rep( "CircleToRectangle", numberOfPoints ), rep( "Rectangle", numberOfPoints ) )
> 
> pointPlotDf <- data.frame( X = allPointsX, Y = allPointsY, Type = pointsType )
> pointPlot <- ggplot( data = pointPlotDf ) +
>              geom_point( aes( x = X, y = Y, colour = Type, shape = Type ) ) +
>              geom_path( aes( x = X, y = Y, colour = Type ), alpha = 0.25 ) +
>              ggtitle( "Warping between endpoints" )
> 
> pointPlot
> 
> # Warp the endpoints (i.e., rectangle and circle) to the middle (i.e., square)
> 
> forwardXfrm <- createAntsrTransform( type = "DisplacementFieldTransform",
>         displacement.field = integrateVelocityField( tv$velocityField, 0, 0.5 ) )
> warpedRectangleToSquare <- applyAntsrTransformToPoint( forwardXfrm, rectangle )
> inverseXfrm <- createAntsrTransform( type = "DisplacementFieldTransform",
>         displacement.field = integrateVelocityField( tv$velocityField, 1.0, 0.5 ) )
> warpedCircleToSquare <- applyAntsrTransformToPoint( inverseXfrm, circle )
> 
> allPointsX <- c( rectangle[,1], square[,1], circle[,1], warpedRectangleToSquare[,1], warpedCircleToSquare[,1]  )
> allPointsY <- c( rectangle[,2], square[,2], circle[,2], warpedRectangleToSquare[,2], warpedCircleToSquare[,2] )
> pointsType <- c( rep( "Rectangle", numberOfPoints ), rep( "Square", numberOfPoints ), rep( "Circle", numberOfPoints ),
>                  rep( "RectangleToSquare", numberOfPoints ), rep( "CircleToSquare", numberOfPoints ) )
> 
> pointPlotDf <- data.frame( X = allPointsX, Y = allPointsY, Type = pointsType )
> pointPlot <- ggplot( data = pointPlotDf ) +
>              geom_point( aes( x = X, y = Y, colour = Type, shape = Type ) ) +
>              geom_path( aes( x = X, y = Y, colour = Type ), alpha = 0.25 ) +
>              ggtitle( "Warping to middle" )
> 
> pointPlot

ANTsR / ANTsPy

Miscellaneous utilities

Bias correction and denoising

Python

>>> import ants
>>>
>>> r16 = ants.image_read(ants.get_ants_data("r16"))
>>> r16_mask = ants.get_mask(r16)
>>> ants.plot(r16, overlay=r16_mask, overlay_alpha = 0.5)
>>> 
>>> # Bias correction
>>> r16_n3 = ants.n3_bias_field_correction2(r16, r16_mask, verbose=True)
>>> r16_n4 = ants.n4_bias_field_correction(r16, r16_mask, verbose=True)
>>> 
>>> # Denoising
>>> r16_denoised = ants.denoise_image(r16_n3, r16_mask, v=1)

R

> library( ANTsR )
>
> r16 <- antsImageRead( getANTsRData( "r16" ) )
> r16Mask <- getMask( r16 )
> plot( r16, r16Mask, alpha = 0.5 )
> 
> # Bias correction
> r16N3 <- n3BiasFieldCorrection2( r16, r16Mask, verbose = TRUE )
> r16N4 <- n4BiasFieldCorrection( r16, r16Mask, verbose = TRUE )
> 
> # Denoising
> r16Denoised <- denoiseImage( r16N3, r16Mask, verbose = TRUE )

Atropos segmentation and quantification

Python

>>> import ants
>>>
>>> r16 = ants.image_read(ants.get_ants_data("r16"))
>>> r16_mask = ants.get_mask(r16)
>>> 
>>> # Atropos
>>> r16_atropos = ants.atropos(r16, r16_mask)
>>> ants.plot(r16, overlay=r16_atropos['segmentation'], overlay_alpha=0.5)
>>> 
>>> # Label measures
>>> ants.label_geometry_measures(r16_atropos['segmentation'])
   Label  VolumeInMillimeters  SurfaceAreaInMillimetersSquared  Eccentricity  ...  BoundingBoxUpper_y  IntegratedIntensity  WeightedCentroid_x  WeightedCentroid_y
0      1               2328.0                      2018.412135      0.579791  ...                 212                 2328          132.506014          131.398196
1      2               7319.0                      3645.703919      0.540221  ...                 211                14638          124.857084          125.767045
2      3               8370.0                      2155.330184      0.646684  ...                 209                25110          124.823775          131.760096

[3 rows x 17 columns]

R

> library( ANTsR )
>
> r16 <- antsImageRead( getANTsRData( "r16" ) )
> r16Mask <- getMask( r16 )
> 
> # Atropos
> r16Atropos <- atropos( r16, r16Mask )
> plot( r16, r16Mask, alpha = 0.5 )
> 
> # Label measures
> labelGeometryMeasures( r16Atropos$segmentation )
  Label VolumeInMillimeters SurfaceAreaInMillimetersSquared Eccentricity
1     1                2331                        2020.078    0.5797638
2     2                7311                        3646.259    0.5395833
3     3                8375                        2155.886    0.6470552
  Elongation Orientation Centroid_x Centroid_y AxesLength_x AxesLength_y
1   1.227318    1.565385   132.5324   131.4033     137.0855     168.2476
2   1.187744    1.497297   124.8014   125.7078     148.7605     176.6894
3   1.311573    1.447841   124.8623   131.8048     123.3535     161.7871
  BoundingBoxLower_x BoundingBoxUpper_x BoundingBoxLower_y BoundingBoxUpper_y
1                 57                195                 44                212
2                 58                194                 45                211
3                 60                192                 49                209
  IntegratedIntensity WeightedCentroid_x WeightedCentroid_y
1                2331           132.5324           131.4033
2               14622           124.8014           125.7078
3               25125           124.8623           131.8048

Histogram matching

The following examples demonstrate a variant of the well-known Nyul histogram standardization algorithm that permits masking and uses a B-spline approximation technique for histogram fitting. For the ITK implementation, use histogram_matching(...) in ANTsPy or histogramMatching(...) in ANTsR.

Python

>>> import ants
>>> from matplotlib import pyplot
>>> 
>>> source_image = ants.image_read(ants.get_ants_data("r16"))
>>> source_mask = ants.get_mask(source_image)
>>> 
>>> reference_image = ants.image_read(ants.get_ants_data("r64"))
>>> reference_mask = ants.get_mask(reference_image)
>>> 
>>> tx_source_image = ants.histogram_match_image2(source_image, reference_image, 
...                                               source_mask, reference_mask,
...                                               match_points=64,
...                                               transform_domain_size=255)
>>>
>>> source = source_image[source_mask != 0].flatten()
>>> reference = reference_image[reference_mask != 0].flatten()
>>> tx_source = tx_source_image[source_mask != 0].flatten()
>>>
>>> pyplot.hist([source, reference, tx_source], 64, alpha=0.5, label=["Source", "Ref", "Tx"])
>>> pyplot.legend(loc="upper right")
>>> pyplot.show()

R

> library( ANTsR )
> library( ggplot2 )
> 
> sourceImage <- antsImageRead( getANTsRData( "r16" ) )
> referenceImage <- antsImageRead( getANTsRData( "r64" ) )
> sourceMask <- getMask( sourceImage )
> referenceMask <- getMask( referenceImage )
> 
> txSourceImage <- histogramMatchImage2( txSourceImage, referenceImage, 
>                                        sourceMask, referenceMask,
>                                        matchPoints = 64, 
>                                        transformDomainSize = 255 )
> 
> sourceHist <- as.vector( sourceImage[sourceMask != 0] )
> referenceHist <- as.vector( referenceImage[referenceMask != 0] )
> txSourceHist <- as.vector( txSourceImage[sourceMask != 0] )
>
> histData <- data.frame( values = c( sourceHist, referenceHist, txSourceHist ),
>                         imageType = c( rep( "Src", length( sourceHist ) ),  
>                                        rep( "Ref", length( referenceHist ) ),   
>                                        rep( "Tx", length( txSourceHist ) ) ) 
>                       )
>
> histPlot <- ggplot( histData, aes( x = values, color = imageType ) ) +
>                     geom_density( alpha = 0.5, position = "identity")
> histPlot

ANTsR / ANTsPy

Deformation-based morphometry

Python

>>> import ants
>>> import numpy as np
>>> import pandas as pd
>>> import random
>>> import statsmodels
>>> 
>>> # Build template
>>> r16 = ants.image_read(ants.get_ants_data('r16'))
>>> r27 = ants.image_read(ants.get_ants_data('r27'))
>>> r30 = ants.image_read(ants.get_ants_data('r30'))
>>> r62 = ants.image_read(ants.get_ants_data('r62'))
>>> r64 = ants.image_read(ants.get_ants_data('r64'))
>>> r85 = ants.image_read(ants.get_ants_data('r85'))
>>> 
>>> rimage_list = [r16, r27, r30, r62, r64, r85]
>>> number_of_subjects = len(rimage_list)
>>> rtemplate = ants.build_template(initial_template=None, image_list=rimage_list, type_of_transform="antsRegistrationSyNQuick[s]")
>>> 
>>> rtemplate_seg = ants.kmeans_segmentation(rtemplate, 3)['segmentation']
>>> rtemplate_gm = ants.threshold_image(rtemplate_seg, 2, 2, 1, 0)
>>> 
>>> # Generate log jacobian images and get brain volumes for later use
>>> log_jacobian_image_list = list()
>>> brain_volume = list()
>>> for i in range(number_of_subjects):
>>>     print("Registration subject", i, "to template.")
>>>     reg = ants.registration(fixed=rtemplate, moving=rimage_list[i], type_of_transform="antsRegistrationSyNQuick[b]")
>>>     log_jacobian_image_list.append(ants.create_jacobian_determinant_image(rtemplate, reg['fwdtransforms'][0], do_log=True))     
>>>     geom = ants.label_geometry_measures(ants.get_mask(rimage_list[i]))
>>>     brain_volume.append(geom['VolumeInMillimeters'][0])
>>> 
>>> # Create matrix to facilitate statistical analysis
>>> log_jacobian = ants.image_list_to_matrix(log_jacobian_image_list, rtemplate_gm)
>>> 
>>> # Create simulated outcomes and covariates
>>> neuro_condition = np.zeros((number_of_subjects,))
>>> age = np.zeros((number_of_subjects,))
>>> sex = np.zeros((number_of_subjects,))
>>> for i in range(number_of_subjects):
>>>     neuro_condition[i] = random.random()
>>>     age[i] = random.randint(40, 70)
>>>     sex[i] = random.randint(0, 1)
>>> data = {'neuro_condition' : neuro_condition, 'age' : age, 'sex' : sex, 'brain_volume' : brain_volume}
>>> covariates = pd.DataFrame(data)
>>> rformula = "neuro_condition ~ age + sex + brain_volume + log_jacobian"
>>> 
>>> # Run image-based regression
>>> dbm = ants.ilr(covariates, {"log_jacobian" : log_jacobian}, rformula)
>>> 
>>> # FDR correction
>>> log_jacobian_p_values = dbm['pValues']['pval_log_jacobian']
>>> log_jacobian_q_values = statsmodels.stats.multitest.fdrcorrection(log_jacobian_p_values, alpha=0.05, method='poscorr', is_sorted=False)[1]
>>> 
>>> # Create overlay q-value image
>>> log_jacobian_q_values_image = ants.matrix_to_images(np.reshape(log_jacobian_q_values, (1, len(log_jacobian_q_values))), rtemplate_gm)

R

> library( ANTsR )
> 
> # Build template
> r16 <- antsImageRead( getANTsRData( 'r16' ) )
> r27 <- antsImageRead( getANTsRData( 'r27' ) )
> r30 <- antsImageRead( getANTsRData( 'r30' ) )
> r62 <- antsImageRead( getANTsRData( 'r62' ) )
> r64 <- antsImageRead( getANTsRData( 'r64' ) )
> r85 <- antsImageRead( getANTsRData( 'r85' ) )
> 
> rimageList <- list( r16, r27, r30, r62, r64, r85 )
> numberOfSubjects <- length( rimageList )
> initialAverage <- antsAverageImages( rimageList )
> rtemplate <- buildTemplate( initialTemplate = initialAverage, imgList = rimageList, typeofTransform = "antsRegistrationSyNQuick[s]", verbose = TRUE )
> 
> rtemplateSeg <- thresholdImage( rtemplate, "Kmeans", 3 )
> rtemplateGM <- thresholdImage( rtemplateSeg, 2, 2, 1, 0 )
> 
> # Generate log jacobian images and get brain volumes for later use
> logJacobianImageList <- list()
> brainVolume <- c()
> for( i in seq.int( numberOfSubjects ) )
>   {
>   cat( "Registration subject", i, "to template.\n")
>   reg <- antsRegistration( fixed = rtemplate, moving = rimageList[[i]], typeofTransform = "antsRegistrationSyNQuick[b]" )
>   logJacobianImageList[[i]] <- createJacobianDeterminantImage( rtemplate, reg$fwdtransforms[[1]], doLog = TRUE )
> 
>   geom <- labelGeometryMeasures( getMask( rimageList[[i]] ) )
>   brainVolume[i] <- geom$VolumeInMillimeters[1] 
>   }
> 
> # Create matrix to facilitate statistical analysis
> logJacobian <- imageListToMatrix( logJacobianImageList, rtemplateGM )
> 
> # Create simulated outcomes and covariates
> neuroCondition <- rnorm( numberOfSubjects )
> age <- sample( 40:70, numberOfSubjects, replace = TRUE )
> sex <- sample( c( 0, 1 ), numberOfSubjects, replace = TRUE )
> covariates <- data.frame( neuroCondition = neuroCondition, age = age, sex = sex, brainVolume = brainVolume)
> rformula <- "neuroCondition ~ age + sex + brainVolume + logJacobian"
> 
> # Run image-based regression
> dbm <- ilr( covariates, list( logJacobian = logJacobian ), rformula )
>
# FDR correction
> logJacobianPvalues <- dbm$pValue[5,]
> logJacobianQvalues <- p.adjust( logJacobianPvalues, method = "fdr" )
> 
> # Create overlay q-value image
> logJacobianQvaluesImage <- matrixToImages( matrix( data = logJacobianQvalues, nrow = 1 ), rtemplateGM )

ANTsXNet

Installation

  • Python: git clone https://github.com/ANTsX/ANTsPyNet; cd ANTsPyNet; python3 setup.py install
  • R: devtools::install_github( "ANTsX/ANTsRNet" )

A couple important notes

  • Primarily geared towards applications
  • No need for a GPU

ANTsXNet

Brain applications

  • T1 (multiple "flavors")
  • T2
  • T2 star
  • T1/T2 infant
  • mean bold
  • FA

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> 
>>> # ANTs-flavored
>>> seg = antspynet.brain_extraction(t1, modality="t1", verbose=True)
>>> ants.plot(t1, overlay=seg, overlay_alpha=0.5)
>>>
>>> # FreeSurfer-flavored
>>> seg = antspynet.brain_extraction(t1, modality="t1nobrainer", verbose=True)
>>> ants.plot(t1, overlay=seg, overlay_alpha=0.5)
>>>
>>> # Combined
>>> seg = antspynet.brain_extraction(t1, modality="t1combined", verbose=True)
>>> ants.plot(t1, overlay=seg, overlay_alpha=0.5)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> 
> # ANTs-flavored
> seg <- brainExtraction( t1, modality = "t1", verbose = TRUE )
> plot( t1, seg, alpha = 0.5 )
>
> # FreeSurfer-flavored
> seg <- brainExtraction( t1, modality = "t1nobrainer", verbose = TRUE )
> plot( t1, seg, alpha = 0.5 )
>
> # Combined
> seg <- brainExtraction( t1, modality = "t1combined", verbose = TRUE )
> plot( t1, seg, alpha = 0.5 )

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> seg = antspynet.deep_atropos(t1, verbose=True)
>>> ants.plot(t1, overlay=seg['segmentation_image'], overlay_alpha=0.75)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> seg <- deepAtropos( t1, verbose = TRUE )
> plot( t1, seg$segmentationImage, alpha = 0.75 )

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> kk = antspynet.cortical_thickness(t1, verbose=True)
>>> ants.plot(t1, overlay=kk['thickness_image'], overlay_alpha=0.75)
>>>
>>> # Also see antspynet.longitudinal_cortical_thickness(...)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> kk <- corticalThickness( t1, verbose = TRUE )
> plot( t1, kk$thicknessImage, alpha = 0.75 )
>
> # Also see longitudinalCorticalThickness(...)

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> dkt = antspynet.desikan_killiany_tourville_labeling(t1, do_lobar_parcellation=True, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> dkt <- desikanKillianyTourvilleLabeling( t1, doLobarParcellation = TRUE, verbose = TRUE )
  • T1
  • T1/T2

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> df = antspynet.deep_flash(t1, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> df <- deepFlash( t1, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> hipp = antspynet.hippmapp3r_segmentation(t1, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> hipp <- hippMapp3rSegmentation( t1, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> age = antspynet.brain_age(t1, number_of_simulations=3, sd_affine=0.01, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> age <- brainAge( t1, numberOfSimulations = 3, sdAffine = 0.01, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> seg = antspynet.claustrum_segmentation(t1, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> seg <- claustrumSegmentation( t1, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('kirby'))
>>> seg = antspynet.hypothalamus_segmentation(t1, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'kirby' ) )
> seg <- hypothalamusSegmentation( t1, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> t1_file = tf.keras.utils.get_file(fname="t1.nii.gz", origin="https://figshare.com/ndownloader/files/40251796", force_download=True)
>>> t1 = ants.image_read(t1_file)
>>> flair_file = tf.keras.utils.get_file(fname="flair.nii.gz", origin="https://figshare.com/ndownloader/files/40251793", force_download=True)
>>> flair = ants.image_read(flair_file)
>>>
>>> wmh = antspynet.sysu_media_wmh_segmentation(flair, t1, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1File <- tensorflow::tf$keras$utils$get_file( fname = "t1.nii.gz", origin = "https://figshare.com/ndownloader/files/40251796", force_download = TRUE )
> t1 <- antsImageRead( t1File )
> flairFile <- tensorflow::tf$keras$utils$get_file( fname = "flair.nii.gz", origin = "https://figshare.com/ndownloader/files/40251793", force_download = TRUE )
> flair <- antsImageRead( flairFile )
> 
> wmh <- sysuMediaWmhSegmentation( flair, t1, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> t1_file = tf.keras.utils.get_file(fname="t1.nii.gz", origin="https://figshare.com/ndownloader/files/40251796", force_download=True)
>>> t1 = ants.image_read(t1_file)
>>> flair_file = tf.keras.utils.get_file(fname="flair.nii.gz", origin="https://figshare.com/ndownloader/files/40251793", force_download=True)
>>> flair = ants.image_read(flair_file)
>>>
>>> wmh = antspynet.hypermapp3r_segmentation(t1, flair, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1File <- tensorflow::tf$keras$utils$get_file( fname = "t1.nii.gz", origin = "https://figshare.com/ndownloader/files/40251796", force_download = TRUE )
> t1 <- antsImageRead( t1File )
> flairFile <- tensorflow::tf$keras$utils$get_file( fname = "flair.nii.gz", origin = "https://figshare.com/ndownloader/files/40251793", force_download = TRUE )
> flair <- antsImageRead( flairFile )
> 
> wmh <- hyperMapp3rSegmentation( t1, flair, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> t1_file = tf.keras.utils.get_file(fname="t1.nii.gz", origin="https://figshare.com/ndownloader/files/40251796", force_download=True)
>>> t1 = ants.image_read(t1_file)
>>> flair_file = tf.keras.utils.get_file(fname="flair.nii.gz", origin="https://figshare.com/ndownloader/files/40251793", force_download=True)
>>> flair = ants.image_read(flair_file)
>>>
>>> wmh = antspynet.shiva_wmh_segmentation(flair, t1, which_model = "all", verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1File <- tensorflow::tf$keras$utils$get_file( fname = "t1.nii.gz", origin = "https://figshare.com/ndownloader/files/40251796", force_download = TRUE )
> t1 <- antsImageRead( t1File )
> flairFile <- tensorflow::tf$keras$utils$get_file( fname = "flair.nii.gz", origin = "https://figshare.com/ndownloader/files/40251793", force_download = TRUE )
> flair <- antsImageRead( flairFile )
> 
> wmh <- shivaWmhSegmentation( flair, t1, whichModel = "all", verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> t1_file = tf.keras.utils.get_file(fname="t1.nii.gz", origin="https://figshare.com/ndownloader/files/40251796", force_download=True)
>>> t1 = ants.image_read(t1_file)
>>> t1 = ants.resample_image(t1, (240, 240, 64), use_voxels=True)
>>> flair_file = tf.keras.utils.get_file(fname="flair.nii.gz", origin="https://figshare.com/ndownloader/files/40251793", force_download=True)
>>> flair = ants.image_read(flair_file)
>>> flair = ants.resample_image(flair, (240, 240, 64), use_voxels=True)
>>>
>>> wmh = antspynet.wmh_segmentation(flair, t1, use_combined_model=True, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1File <- tensorflow::tf$keras$utils$get_file( fname = "t1.nii.gz", origin = "https://figshare.com/ndownloader/files/40251796", force_download = TRUE )
> t1 <- antsImageRead( t1File )
> t1 <- resampleImage( t1, c( 240, 240, 64 ), useVoxels = TRUE )
> flairFile <- tensorflow::tf$keras$utils$get_file( fname = "flair.nii.gz", origin = "https://figshare.com/ndownloader/files/40251793", force_download = TRUE )
> flair <- antsImageRead( flairFile )
> flair <- resampleImage( flair, c( 240, 240, 64 ), useVoxels = TRUE )
> 
> wmh <- wmhSegmentation( t1, flair, useCombinedModel = TRUE, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> t1_file = tf.keras.utils.get_file(fname="pvs_t1.nii.gz", origin="https://figshare.com/ndownloader/files/48675367", force_download=True)
>>> t1 = ants.image_read(t1_file)
>>> flair_file = tf.keras.utils.get_file(fname="pvs_flair.nii.gz", origin="https://figshare.com/ndownloader/files/48675352", force_download=True)
>>> flair = ants.image_read(flair_file)
>>>
>>> pvs = antspynet.shiva_pvs_segmentation(t1, flair, which_model = "all", verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1File <- tensorflow::tf$keras$utils$get_file( fname = "pvs_t1.nii.gz", origin = "https://figshare.com/ndownloader/files/48675367", force_download = TRUE )
> t1 <- antsImageRead( t1File )
> flairFile <- tensorflow::tf$keras$utils$get_file( fname = "pvs_flair.nii.gz", origin = "https://figshare.com/ndownloader/files/48675352", force_download = TRUE )
> flair <- antsImageRead( flairFile )
> 
> pvs <- shivaPvsSegmentation( t1, flair, whichModel = "all", verbose = TRUE )

Cerebellum morphology

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data("mprage_hippmapp3r"))
>>>
>>> # Computing the thickness image will take relatively more time than the other	
>>> # outputs so the recommendation would be to run it initially without computing the	
>>> # thickness image to just get a sense of the application.	
>>>
>>> cereb = antspynet.cerebellum_morphology(t1, compute_thickness_image=True, verbose=True)	
>>> 	
>>> # possible refinement with cerebellum estimate (comment out since it's not needed for this image).	
>>> # mask = ants.threshold_image(cereb['cerebellum_probability_image'], 0.5, 1, 1, 0)	
>>> # cereb = antspynet.cerebellum_morphology(t1, cerebellum_mask=mask, verbose=True)	
>>> 
>>> ants.image_write(cereb['cerebellum_probability_image'], "cerebellum_probability_mask.nii.gz")	
>>> ants.image_write(cereb['thickness_image'], "kk.nii.gz")	
>>> ants.image_write(cereb['parcellation_segmentation_image'], "parcelation.nii.gz")	

R

> library( ANTsR )	
> library( ANTsRNet )	
>	
> t1 <- antsImageRead( getANTsXNetData( "mprage_hippmapp3r" ) )	
>	
> # Computing the thickness image will take relatively more time than the other	
> # outputs so the recommendation would be to run it initially without computing the	
> # thickness image to just get a sense of the application.	
> cereb <- cerebellumMorphology( t1, computeThicknessImage = TRUE, verbose = TRUE )	
>	
> # possible refinement with cerebellum estimate (comment out since it's not needed for this image).	
> # mask <- thresholdImage( cereb$cerebellumProbabilityImage, 0.5, 1, 1, 0 )	
> # cereb <- cerebellumMorphology( t1, cerebellumMask = mask, verbose = TRUE )	
> 	
> # Write output to disk	
> antsImageWrite( cereb$cerebellumProbabilityImage, "cerebellumProbabilityMask.nii.gz" )	
> antsImageWrite( cereb$tissueSegmentationImage, "tissue.nii.gz" )	
> antsImageWrite( cereb$thicknessImage, "kk.nii.gz" )	
> antsImageWrite( cereb$parcellationSegmentationImage, "parcelation.nii.gz" )	

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> flair_file = tf.keras.utils.get_file(fname="flair.nii.gz", origin="https://figshare.com/ndownloader/files/42385077", force_download=True)
>>> flair = ants.image_read(flair_file)
>>> t1_file = tf.keras.utils.get_file(fname="t1.nii.gz", origin="https://figshare.com/ndownloader/files/42385071", force_download=True)
>>> t1 = ants.image_read(t1_file)
>>> t1_contrast_file = tf.keras.utils.get_file(fname="t1_contrast.nii.gz", origin="https://figshare.com/ndownloader/files/42385068", force_download=True)
>>> t1_contrast = ants.image_read(t1_contrast_file)
>>> t2_file = tf.keras.utils.get_file(fname="t2.nii.gz", origin="https://figshare.com/ndownloader/files/42385074", force_download=True)
>>> t2 = ants.image_read(t2_file)
>>>
>>> bt = antspynet.brain_tumor_segmentation(flair, t1, t1_contrast, t2, patch_stride_length=32, verbose=True)
>>> # ants.image_write(bt['segmentation_image'], "brain_tumor_segmentation.nii.gz")

R

> library( ANTsR )
> library( ANTsRNet )
>
> flairFile <- tensorflow::tf$keras$utils$get_file( fname = "flair.nii.gz", origin = "https://figshare.com/ndownloader/files/42385077", force_download = TRUE )
> flair <- antsImageRead( flairFile )
> t1File <- tensorflow::tf$keras$utils$get_file( fname = "t1.nii.gz", origin = "https://figshare.com/ndownloader/files/42385071", force_download = TRUE )
> t1 <- antsImageRead( t1File )
> t1ContrastFile <- tensorflow::tf$keras$utils$get_file( fname = "t1_contrast.nii.gz", origin = "https://figshare.com/ndownloader/files/42385068", force_download = TRUE )
> t1Contrast <- antsImageRead( t1ContrastFile )
> t2File <- tensorflow::tf$keras$utils$get_file( fname = "t2.nii.gz", origin = "https://figshare.com/ndownloader/files/42385074", force_download = TRUE )
> t2 <- antsImageRead( t2File )
> 
> bt <- brainTumorSegmentation( flair, t1, t1Contrast, t2, patchStrideLength = 32, verbose = TRUE )
> #antsImageWrite( bt$segmentationImage, "brainTumorSegmentation.nii.gz" )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> mmbop_file = tf.keras.utils.get_file(fname="mra.nii.gz", origin="https://figshare.com/ndownloader/files/46406755", force_download=True)
>>> mra = ants.image_read(mra_file)
>>> vessels = antspynet.mra_brain_vessel_segmentation(mra, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> mmbopFile <- tensorflow::tf$keras$utils$get_file( fname = "mra.nii.gz", origin = "https://figshare.com/ndownloader/files/46406755", force_download = TRUE )
> mra <- antsImageRead( mmbopFile )
> vessels <- mraBrainVesselSegmentation( mra, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> t1_file = tf.keras.utils.get_file(fname="t1w_with_lesion.nii.gz", origin="https://figshare.com/ndownloader/files/44053868", force_download=True)
>>> t1 = ants.image_read(t1_file)
>>>
>>> probability_mask = antspynet.lesion_segmentation(t1, do_preprocessing=True, verbose=True)
>>> ants.image_write(probability_mask, "lesion_probability_mask.nii.gz")

R

> library( ANTsR )
> library( ANTsRNet )
> library( tensorflow )
>
> t1File <- tensorflow::tf$keras$utils$get_file( fname = "t1w_with_lesion.nii.gz", origin = "https://figshare.com/ndownloader/files/44053868", force_download = TRUE )
> t1 <- antsImageRead( t1File )
>
> probabilityMask <- lesionSegmentation( t1, doPreprocessing = TRUE, verbose = TRUE )
> antsImageWrite( probabilityMask, "lesion_probability_mask.nii.gz" )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> t1_file = tf.keras.utils.get_file(fname="t1w_with_lesion.nii.gz", origin="https://figshare.com/ndownloader/files/44053868", force_download=True)
>>> t1 = ants.image_read(t1_file)
>>>
>>> probability_mask = antspynet.lesion_segmentation(t1, do_preprocessing=True, verbose=True)
>>> lesion_mask = ants.threshold_image(probability_mask, 0.5, 1.1, 1, 0)
>>> t1_inpainted = antspynet.whole_head_inpainting(t1, roi_mask=lesion_mask, modality="t1", mode="axial", verbose=True)
>>> ants.image_write(t1_inpainted, "t1_repaired.nii.gz")

R

> library( ANTsR )
> library( ANTsRNet )
> library( tensorflow )
>
> t1File <- tensorflow::tf$keras$utils$get_file( fname = "t1w_with_lesion.nii.gz", origin = "https://figshare.com/ndownloader/files/44053868", force_download = TRUE )
> t1 <- antsImageRead( t1File )
>
> probabilityMask <- lesionSegmentation( t1, doPreprocessing = TRUE, verbose = TRUE )
> lesionMask <- thresholdImage( probabilityMask, 0.5, 1.1, 1, 0 )
> t1Inpainted <- wholeHeadInpainting( t1, roiMask = lesionMask, modality = "t1", mode = "axial", verbose = TRUE )
> antsImageWrite( t1Inpainted, "t1_repaired.nii.gz" )

Lung applications

  • CT (left, right, and airways)
  • Proton (left, right)
  • Ventilation-based whole lung segmentation
  • X-ray (left, right)

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> # CT
>>> ct_file = tf.keras.utils.get_file(fname="ctLung.nii.gz", origin="https://figshare.com/ndownloader/files/42934234", force_download=True)
>>> ct = ants.image_read(ct_file)
>>> lung_ex = antspynet.lung_extraction(ct, modality="ct", verbose=True)
>>> ants.plot(ants.iMath_normalize(ct), lung_ex['segmentation_image'], alpha=0.5)
>>>
>>> # Proton (whole lung)
>>> proton_file = tf.keras.utils.get_file(fname="protonLung.nii.gz", origin="https://figshare.com/ndownloader/files/42934228", force_download=True)
>>> proton = ants.image_read(proton_file)
>>> lung_ex = antspynet.lung_extraction(proton, modality="proton", verbose=True)
>>> ants.plot(ants.iMath_normalize(proton), lung_ex['segmentation_image'], alpha=0.5)
>>>
>>> # Proton (lobes)
>>> lung_ex = antspynet.lung_extraction(proton, modality="protonLobes", verbose=True)
>>> ants.plot(ants.iMath_normalize(proton), lung_ex['segmentation_image'], alpha=0.5)
>>>
>>> # Lobes from lung mask
>>> lung_mask = ants.threshold_image(lung_ex['segmentation_image'], 0, 0, 0, 1 )
>>> lung_ex = antspynet.lung_extraction(lung_mask, modality="maskLobes", verbose=True)
>>> ants.plot(ants.iMath_normalize(proton), lung_ex['segmentation_image'], alpha=0.5)
>>>
>>> # Chest x-ray
>>> cxr_file = tf.keras.utils.get_file(fname="cxr.nii.gz", origin="https://figshare.com/ndownloader/files/42934237", force_download=True)
>>> cxr = ants.image_read(cxr_file)
>>> lung_ex = antspynet.lung_extraction(cxr, modality="xray", verbose=True)
>>> ants.plot(ants.iMath_normalize(cxr), lung_ex['segmentation_image'], alpha=0.5)
>>>
>>> # Ventilation
>>> ventilation_file = tf.keras.utils.get_file(fname="ventilationLung.nii.gz", origin="https://figshare.com/ndownloader/files/42934231", force_download=True)
>>> ventilation = ants.image_read(ventilation_file)
>>> lung_ex = antspynet.lung_extraction(ventilation, modality="ventilation", verbose=True)
>>> ants.plot(ants.iMath_normalize(ventilation), ants.threshold_image(lung_ex, 0.5, 1, 1, 0), alpha=0.5)

R

> library( ANTsR )
> library( ANTsRNet )
>
> # CT
> ctFile <- tensorflow::tf$keras$utils$get_file( fname = "ctLung.nii.gz", origin = "https://figshare.com/ndownloader/files/42934234", force_download = TRUE )
> ct <- antsImageRead( ctFile )
> lungEx <- lungExtraction( ct, modality = "ct", verbose = TRUE )
> plot( ct, lungEx$segmentationImage, alpha = 0.5 )
>
> # Proton (whole lung)
> protonFile <- tensorflow::tf$keras$utils$get_file( fname = "protonLung.nii.gz", origin = "https://figshare.com/ndownloader/files/42934228", force_download = TRUE )
> proton <- antsImageRead( protonFile )
> lungEx <- lungExtraction( proton, modality = "proton", verbose = TRUE )
> plot( proton, lungEx$segmentationImage, alpha = 0.5 )
>
> # Proton (lobes)
> lungEx <- lungExtraction( proton, modality = "protonLobes", verbose = TRUE )
> plot( proton, lungEx$segmentationImage, alpha = 0.5 )
>
> # Lobes from lung mask
> lungMask <- thresholdImage( lungEx$segmentationImage, 0, 0, 0, 1 )
> lungEx <- lungExtraction( lungMask, modality = "maskLobes", verbose = TRUE )
> plot( proton, lungEx$segmentationImage, alpha = 0.5 )
>
> # Chest x-ray
> cxrFile <- tensorflow::tf$keras$utils$get_file( fname = "cxr.nii.gz", origin = "https://figshare.com/ndownloader/files/42934237", force_download = TRUE )
> cxr <- antsImageRead( cxrFile )
> lungEx <- lungExtraction( cxr, modality = "xray", verbose = TRUE )
> plot( cxr, lungEx$segmentationImage, alpha = 0.5 )
>
> # Ventilation
> ventilationFile <- tensorflow::tf$keras$utils$get_file( fname = "ventilationLung.nii.gz", origin = "https://figshare.com/ndownloader/files/42934231", force_download = TRUE )
> ventilation <- antsImageRead( ventilationFile )
> lungEx <- lungExtraction( ventilation, modality = "ventilation", verbose = TRUE )
> plot( ventilation, thresholdImage( lungEx, 0.5, 1, 1, 0 ), alpha = 0.5 )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> # Process the proton image to get the lung mask
>>> proton_file = tf.keras.utils.get_file(fname="protonLung.nii.gz", origin="https://figshare.com/ndownloader/files/42934228", force_download=True)
>>> proton = ants.image_read(proton_file)
>>> lung_ex = antspynet.lung_extraction(proton, modality="proton", verbose=True)
>>> lung_mask = ants.threshold_image(lung_ex['segmentation_image'], 0, 0, 0, 1 )
>>>
>>> # Use 'el bicho' to label the four ventilation levels
>>> ventilation_file = tf.keras.utils.get_file(fname="ventilationLung.nii.gz", origin="https://figshare.com/ndownloader/files/42934231", force_download=True)
>>> ventilation = ants.image_read(ventilation_file)
>>> eb = antspynet.el_bicho(ventilation, lung_mask, verbose=True)
>>> ants.plot(ants.iMath_normalize(ventilation), eb['segmentation_image'], alpha=0.5)

R

> library( ANTsR )
> library( ANTsRNet )
>
> # Process the proton image to get the lung mask
> protonFile <- tensorflow::tf$keras$utils$get_file( fname = "protonLung.nii.gz", origin = "https://figshare.com/ndownloader/files/42934228", force_download = TRUE )
> proton <- antsImageRead( protonFile )
> lungEx <- lungExtraction( proton, modality = "proton", verbose = TRUE )
> lungMask <- thresholdImage( lungEx$segmentationImage, 0, 0, 0, 1 )
>
> # Use 'el bicho' to label the four ventilation levels
> ventilationFile <- tensorflow::tf$keras$utils$get_file( fname = "ventilationLung.nii.gz", origin = "https://figshare.com/ndownloader/files/42934231", force_download = TRUE )
> ventilation <- antsImageRead( ventilationFile )
> eb <- elBicho( ventilation, lungMask, verbose = TRUE )
> plot( ventilation, eb$segmentationImage, alpha = 0.5 )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> ct_file = tf.keras.utils.get_file(fname="ctLung.nii.gz", origin="https://figshare.com/ndownloader/files/42934234", force_download=True)
>>> ct = ants.image_read(ct_file)
>>> arteries = antspynet.lung_pulmonary_artery_segmentation(ct, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> ctFile <- tensorflow::tf$keras$utils$get_file( fname = "ctLung.nii.gz", origin = "https://figshare.com/ndownloader/files/42934234", force_download = TRUE )
> ct <- antsImageRead( ctFile )
> arteries <- lungPulmonaryArterySegmentation( ct, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>> import tensorflow as tf
>>>
>>> ct_file = tf.keras.utils.get_file(fname="ctLung.nii.gz", origin="https://figshare.com/ndownloader/files/42934234", force_download=True)
>>> ct = ants.image_read(ct_file)
>>> airways = antspynet.lung_airway_segmentation(ct, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> ctFile <- tensorflow::tf$keras$utils$get_file( fname = "ctLung.nii.gz", origin = "https://figshare.com/ndownloader/files/42934234", force_download = TRUE )
> ct <- antsImageRead( ctFile )
> airways <- lungAirwaySegmentation( ct, verbose = TRUE )

Python

>>> import ants
>>> import antspynet
>>> import deepsimlr
>>>
>>> # This particular subject is diagnosed as having cardiomegaly 
>>> cxr_file = tf.keras.utils.get_file(fname="cxr.nii.gz", origin="https://figshare.com/ndownloader/files/42429156", force_download=True)
>>> cxr = ants.image_read(cxr_file)
>>> 
>>> # Pytorch version adapted from https://github.com/jrzech/reproduce-chexnet
>>> deepsimlr.chexnet(cxr)
   Atelectasis  Cardiomegaly  Effusion  Infiltration      Mass    Nodule  Pneumonia  Pneumothorax  Consolidation     Edema  Emphysema  Fibrosis  Pleural_Thickening    Hernia
0     0.027259      0.745355  0.042421      0.049225  0.003613  0.026083   0.005946      0.003007       0.005099  0.002063   0.017297  0.059956             0.04354  0.004266
>>>
>>> # Keras implementation in ANTsXNet (original)
>>> antspynet.chexnet(cxr, use_antsxnet_variant=False)
   Atelectasis  Cardiomegaly  Effusion  Infiltration      Mass    Nodule  Pneumonia  Pneumothorax  Consolidation     Edema  Emphysema  Fibrosis  Pleural_Thickening   Hernia
0     0.133284       0.64578  0.085462       0.09486  0.019405  0.018975   0.013067      0.021075       0.028001  0.005851   0.018955    0.0555            0.034219  0.00972
>>>
>>> # Keras implementation in ANTsXNet (variation)
>>> antspynet.chexnet(cxr, use_antsxnet_variant=True, include_tuberculosis_diagnosis=True)
   Atelectasis  Cardiomegaly  Effusion  Infiltration      Mass    Nodule  Pneumonia  Pneumothorax  Consolidation     Edema  Emphysema  Fibrosis  Pleural_Thickening    Hernia  Tuberculosis
0     0.069623      0.898534  0.039589      0.055751  0.005931  0.030309   0.006876      0.013999       0.010621  0.005677   0.048858  0.021668            0.013867  0.008546      0.000285

R

> library( ANTsR )
> library( ANTsRNet )
>
> cxrFile <- tensorflow::tf$keras$utils$get_file( fname = "cxr.nii.gz", origin = "https://figshare.com/ndownloader/files/42429156", force_download = TRUE )
> cxr <- antsImageRead( cxrFile )
> 
> # Keras implementation in ANTsXNet (original)
> chexnet( cxr, useANTsXNetVariant = FALSE )
  Atelectasis Cardiomegaly   Effusion Infiltration       Mass     Nodule
1   0.1065275    0.7038919 0.09204669    0.1006764 0.02943837 0.02722431
   Pneumonia Pneumothorax Consolidation       Edema  Emphysema   Fibrosis
1 0.01293034   0.01917528    0.02944043 0.006727949 0.01533213 0.06293546
  Pleural_Thickening     Hernia
1         0.03341614 0.01022097
> 
> # Keras implementation in ANTsXNet (variation)
> chexnet( cxr, useANTsXNetVariant = TRUE, includeTuberculosisDiagnosis = TRUE )
  Atelectasis Cardiomegaly  Effusion Infiltration       Mass     Nodule
1  0.07443799    0.9660769 0.1168589   0.08752631 0.01058979 0.03241358
   Pneumonia Pneumothorax Consolidation      Edema Emphysema  Fibrosis
1 0.01211402   0.02042338    0.01615181 0.01346725  0.052293 0.0357863
  Pleural_Thickening      Hernia Tuberculosis
1          0.0235217 0.008928899 0.0002852035
>

Evaluation on CheXNet testing data

Compare with original implementation

AUC values

                     Pytorch     Keras  ANTsXNet
Atelectasis         0.832871  0.822537  0.863367
Cardiomegaly        0.946268  0.938183  0.962786
Consolidation       0.816950  0.815728  0.855411
Edema               0.912380  0.915611  0.938205
Effusion            0.893323  0.891513  0.916355
Emphysema           0.906865  0.887093  0.920824
Fibrosis            0.834109  0.827323  0.868360
Hernia              0.882469  0.824022  0.919963
Infiltration        0.770469  0.769772  0.805857
Mass                0.875950  0.834911  0.893959
Nodule              0.788035  0.764380  0.823646
Pleural_Thickening  0.828046  0.814205  0.854000
Pneumonia           0.758576  0.751609  0.800324
Pneumothorax        0.880895  0.870271  0.910933


Mouse applications

Mouse brain extraction

Python

>>> import ants
>>> import antspynet
>>>
>>> mouse_t2_file = tf.keras.utils.get_file(fname="mouse.nii.gz", origin="https://figshare.com/ndownloader/files/45289309", force_download=True)
>>> mouse_t2 = ants.image_read(mouse_t2_file)
>>> mouse_t2_n4 = ants.n4_bias_field_correction(mouse_t2, 
                                                rescale_intensities=True,
                                                shrink_factor=2, 
                                                convergence={'iters': [50, 50, 50, 50], 'tol': 0.0}, 
                                                spline_param=20, verbose=True)
>>> mask = antspynet.mouse_brain_extraction(mouse_t2_n4, modality='t2', verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> mouseT2File <- tensorflow::tf$keras$utils$get_file( fname="mouse.nii.gz",  origin = "https://figshare.com/ndownloader/files/45289309", force_download = TRUE )
> mouseT2 <- antsImageRead( mouseT2File )
> mouseT2N4 <- n4BiasFieldCorrection( mouseT2, 
                                      rescaleIntensities = TRUE,
                                      shrinkFactor = 2, 
                                      convergence = list( iters = c( 50, 50, 50, 50 ), tol = 0.0 ), 
                                      splineParam = 20, verbose = TRUE )
>>> mask <- mouseBrainExtraction( mouseT2N4, modality = 't2', verbose = TRUE )

Mouse brain parcellation

Python

>>> import ants
>>> import antspynet
>>>
>>> mouse_t2_file = tf.keras.utils.get_file(fname="mouse.nii.gz", origin="https://figshare.com/ndownloader/files/45289309", force_download=True)
>>> mouse_t2 = ants.image_read(mouse_t2_file)
>>> mouse_t2_n4 = ants.n4_bias_field_correction(mouse_t2, 
                                                rescale_intensities=True,
                                                shrink_factor=2, 
                                                convergence={'iters': [50, 50, 50, 50], 'tol': 0.0}, 
                                                spline_param=20, verbose=True)
>>> parc_nick = antspynet.mouse_brain_parcellation(mouse_t2_n4, 
                                                   mask=None, 
                                                   which_parcellation="nick",      
                                                   return_isotropic_output=True,  
                                                   verbose=True)
>>> parc_tct = antspynet.mouse_brain_parcellation(mouse_t2_n4, 
                                                  mask=None, 
                                                  which_parcellation="tct",      
                                                  return_isotropic_output=True,  
                                                  verbose=True)                                                      

R

> library( ANTsR )
> library( ANTsRNet )
>
> mouseT2File <- tensorflow::tf$keras$utils$get_file( fname="mouse.nii.gz",  origin = "https://figshare.com/ndownloader/files/45289309", force_download = TRUE )
> mouseT2 <- antsImageRead( mouseT2File )
> mouseT2N4 <- n4BiasFieldCorrection( mouseT2, 
                                      rescaleIntensities = TRUE,
                                      shrinkFactor = 2, 
                                      convergence = list( iters = c( 50, 50, 50, 50 ), tol = 0.0 ), 
                                      splineParam = 20, verbose = TRUE )
> parcNick <- mouseBrainParcellation( mouseT2N4, 
                                      mask = NULL,
                                      whichParcellation = 'nick', 
                                      returnIsotropicOutput = TRUE,
                                      verbose = TRUE )
> parcTct <- mouseBrainParcellation( mouseT2N4, 
                                     mask = NULL,
                                     whichParcellation = 'tct', 
                                     returnIsotropicOutput = TRUE,
                                     verbose = TRUE )                                        

Mouse cortical thickness

Python

>>> import ants
>>> import antspynet
>>>
>>> mouse_t2_file = tf.keras.utils.get_file(fname="mouse.nii.gz", origin="https://figshare.com/ndownloader/files/45289309", force_download=True)
>>> mouse_t2 = ants.image_read(mouse_t2_file)
>>> mouse_t2_n4 = ants.n4_bias_field_correction(mouse_t2, 
                                                rescale_intensities=True,
                                                shrink_factor=2, 
                                                convergence={'iters': [50, 50, 50, 50], 'tol': 0.0}, 
                                                spline_param=20, verbose=True)
>>> kk = antspynet.mouse_cortical_thickness(mouse_t2_n4, 
                                            mask=None, 
                                            return_isotropic_output=True,                                    
                                            verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> mouseT2File <- tensorflow::tf$keras$utils$get_file( fname="mouse.nii.gz",  origin = "https://figshare.com/ndownloader/files/45289309", force_download = TRUE )
> mouseT2 <- antsImageRead( mouseT2File )
> mouseT2N4 <- n4BiasFieldCorrection( mouseT2, 
                                      rescaleIntensities = TRUE,
                                      shrinkFactor = 2, 
                                      convergence = list( iters = c( 50, 50, 50, 50 ), tol = 0.0 ), 
                                      splineParam = 20, verbose = TRUE )
>>> kk <- mouseCorticalThickness( mouseT2N4, 
                                  mask = NULL,
                                  returnIsotropicOutput = TRUE,
                                  verbose = TRUE )

General applications

MRI super resolution

Python

>>> import ants
>>> import antspynet
>>>
>>> # This method is deprecated.  Please see https://github.com/stnava/siq.
>>> 
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> t1_lr = ants.resample_image(t1, (4, 4, 4), use_voxels=False)
>>> ants.plot(t1_lr)
>>> t1_sr = antspynet.mri_super_resolution(t1_lr, verbose=True)
>>> ants.plot(t1_sr)

R

> library( ANTsR )
> library( ANTsRNet )
>
> # This method is deprecated.  Please see https://github.com/stnava/siq.
> 
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> t1Lr <- resampleImage( t1, c( 4, 4, 4 ), useVoxels = FALSE )
> plot( t1Lr )
> t1Sr <- mriSuperResolution( t1Lr, verbose = TRUE )
> plot( t1Sr )

No reference image quality assesment using TID

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> qa = antspynet.tid_neural_image_assessment(t1, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> qa <- tidNeuralImageAssessment( t1, verbose = TRUE )

Full reference image quality assessment

Python

>>> import ants
>>> import antspynet 
>>> 
>>> r16 = ants.image_read(ants.get_data("r16"))
>>> r64 = ants.image_read(ants.get_data("r64"))
>>> psnr_value = antspynet.psnr(r16, r64)
>>> ssim_value = antspynet.ssim(r16, r64)

R

> library( ANTsR )
> library( ANTsRNet ) 
> 
> r16 <- antsImageRead( getANTsRData( "r16" ) )
> r64 <- antsImageRead( getANTsRData( "r64" ) )
> psnrValue <- PSNR( as.array( r16 ), as.array( r64 ) )
> ssimValue <- SSIM( as.array( r16 ), as.array( r64 ) )

Data augmentation

Noise

Python

>>> import ants
>>> import antspynet
>>> import random
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> 
>>> t1 = ants.iMath( t1, "Normalize" )
>>> noise_parameters = (0.0, random.uniform(0, 0.05))
>>> t1 = ants.add_noise_to_image(t1, noise_model="additivegaussian", noise_parameters=noise_parameters)
>>>       

R

> library( ANTsR )
> library( ANTsRNet )
> 
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> 
> t1 <- t1 %>% iMath( "Normalize" )
> noiseParameters <- c( 0.0, runif( 1, 0.0, 0.05 ) )
> t1 <- addNoiseToImage( t1, noiseModel = "additivegaussian", noiseParameters = noiseParameters )

Histogram intensity warping

Python

>>> import ants
>>> import antspynet
>>> import random
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> 
>>> break_points = [0.2, 0.4, 0.6, 0.8]
>>> displacements = list()
>>> for b in range(len(break_points)):
>>>     displacements.append(abs(random.gauss(0, 0.05)))
>>>     if random.sample((True, False), 1)[0]:
>>>         displacements[b] *= -1
>>>
>>> t1 = antspynet.histogram_warp_image_intensities(t1, 
                                                    break_points=break_points, 
                                                    clamp_end_points=(True, False), 
                                                    displacements=displacements)

R

> library( ANTsR )
> library( ANTsRNet )
> 
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> 
> breakpoints <- c( 0.2, 0.4, 0.6, 0.8 )
> displacements <- c()
> for( b in seq.int( length( breakpoints ) ) )
>   {
>   displacements[b] <- abs( runif( 1, 0, 0.05 ) )
>   if( sample( c( TRUE, FALSE ), 1 ) )
>     {
>     displacements[b] = displacements[b] * -1
>     }
>   } 
> 
> t1 <- histogramWarpImageIntensities( t1, 
                                       breakPoints = breakpoints, 
                                       clampEndPoints = c( TRUE, FALSE ), 
                                       displacements = displacements )

Simulate bias field

Python

>>> import ants
>>> import numpy as np
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> log_field = antspynet.simulate_bias_field(t1, 
                                              number_of_points=10, 
                                              sd_bias_field=1.0, 
                                              number_of_fitting_levels=2, 
                                              mesh_size=10)
>>> log_field = log_field.iMath("Normalize")
>>> log_field_array = np.power(np.exp(log_field.numpy()), 4)
>>> t1 = t1 * ants.from_numpy(log_field_array, origin=t1.origin, spacing=t1.spacing, direction=t1.direction)

R

> library( ANTsR )
> 
> t1 <- antsImageRead( getANTsXnetData( "mprage_hippmapp3r" ) )
> logField <- simulateBiasField(image, 
                                numberOfPoints = 10, 
                                sdBiasField = 1.0, 
                                numberOfFittingLevels = 2, 
                                meshSize = 10 ) %>% iMath( "Normalize" )
> logField <- ( exp( logField ) )^4
> t1 <- t1 * logField

Random spatial transformations

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>> segmentation = ants.threshold_image(t1, "Otsu", 3)
>>>
>>> data_augmentation = antspynet.randomly_transform_image_data(t1,
             [[t1]],
             [segmentation],
             number_of_simulations=2,
             transform_type='affineAndDeformation',
             sd_affine=0.01,
             deformation_transform_type="bspline",
             number_of_random_points=1000,
             sd_noise=2.0,
             number_of_fitting_levels=4,
             mesh_size=1,
             sd_smoothing=4.0,
             input_image_interpolator='linear',
             segmentation_image_interpolator='nearestNeighbor')
>>>
>>> simulated_image1 = data_augmentation['simulated_images'][0][0]
>>> simulated_segmentation1 = data_augmentation['simulated_segmentation_images'][0]
>>> simulated_image2 = data_augmentation['simulated_images'][1][0]
>>> simulated_segmentation2 = data_augmentation['simulated_segmentation_images'][1]

R

> library( ANTsR )
> library( ANTsRNet )
> 
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> segmentation <- thresholdImage( t1, "Otsu", 3 )
 Otsu Thresh with 3 thresholds
> 
> dataAugmentation <- randomlyTransformImageData( t1,
              list( list( t1 ) ),
              list( segmentation ),
              numberOfSimulations = 2,
              transformType = 'affineAndDeformation',
              sdAffine = 0.01,
              deformationTransformType = "bspline",
              numberOfRandomPoints = 1000,
              sdNoise = 2.0,
              numberOfFittingLevels = 4,
              meshSize = 1,
              sdSmoothing=4.0,
              inputImageInterpolator = 'linear',
              segmentationImageInterpolator = 'nearestNeighbor' )
> 
> simulatedImage1 <- dataAugmentation$simulatedImages[0][0]
> simulatedSegmentation1 <- dataAugmentation$simulatedSegmentationImages[0]
> simulatedImage2 <- dataAugmentation$simulatedImages[1][0]
> simulatedSegmentation2 <- dataAugmentation$simulatedSegmentationImages[1]

Combined

Python

>>> import antspynet
>>> import ants
>>> 
>>> t1s = [[ants.image_read(antspynet.get_antsxnet_data("kirby"))]]
>>> aug = antspynet.data_augmentation(t1s, segmentation_image_list=None, pointset_list=None, 
>>>                                     number_of_simulations=10, 
>>>                                     reference_image=None, 
>>>                                     transform_type='affineAndDeformation', 
>>>                                     noise_model='additivegaussian', 
>>>                                     noise_parameters=(0.0, 0.05), 
>>>                                     sd_simulated_bias_field=0.05, 
>>>                                     sd_histogram_warping=0.05, 
>>>                                     sd_affine=0.05, 
>>>                                     output_numpy_file_prefix=None, 
>>>                                     verbose=True)
>>> for i in range(len(aug['simulated_images'])):
>>>     print("Writing simulated image ", str(i))
>>>     ants.image_write(aug['simulated_images'][i][0], "simulated_image" + str(i) + ".nii.gz")

R

> library( ANTsR )
> library( ANTsRNet )
> 
> t1s <- list( list( antsImageRead( getANTsXNetData( 'kirby' ) ) ) )
> aug <- dataAugmentation( t1s, segmentationImageList = NULL,
              pointsetList = NULL,
              numberOfSimulations = 10, 
              referenceImage = NULL,
              transformType = 'affineAndDeformation',
              sdAffine = 0.05,
              noiseModel = 'additivegaussian',
              noiseParameters = c( 0.0, 0.05 ),
              sdSimulatedBiasField = 1.0,
              sdHistogramWarping = 0.05,
              outputNumpyFilePrefix = NULL,
              verbose = TRUE )
> for( i in seq.int( length( aug$simulatedImages[[0]] ) )
>   {
>   cat( "Writing simulated image", i, "\n" )
>   antsImageWrite( aug$simulatedImages[[i][[0]]], paste0( "simulatedImages", i, ".nii.gz" ) )
>   }
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment