ANTsX Tutorial

a_background.md

ANTsX ecosystem

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

• Installation

• Basics

  • Read, write, and visualize images
  • Image info
  • Conversion to/from native types
  • Arithmetic operations
  • Other operations

• B-spline fitting

• Image registration

• Label image registration

• Template building

• Motion correction

• Fit transforms to paired points

• Velocity flows across point sets

• Deformation-based morphometry

• Miscellaneous utilities

  • Bias correction and denoising
  • Atropos segmentation and quantification
  • Histogram matching

• fMRI-specific

ANTsXNet

• Brain applications

  • Multi-modal brain extraction
  • Deep Atropos (Six-tissue brain segmentation)
  • Cortical thickness
  • Desikan-Killiany-Tourville parcellation
  • Harvard-Oxford Atlas subcortical parcellation
  • DeepFLASH (medial temporal lobe parcellation)
  • Hippmapp3r (hippocampal segmentation)
  • Brain AGE
  • Claustrum segmentation
  • Hypothalamus segmentation
  • Cerebellum morphology
  • White matter hyperintensities segmentation
    • SYSU
    • Hypermapp3r
    • SHIVA
    • ANTsXNet
  • Perivascular spaces segmentation (SHIVA)
  • Brain tumor segmentation
  • MRA-TOF vessel segmentation
  • Lesion segmentation (WIP)
  • Whole head inpainting

• Lung applications

  • Lung extraction
  • Functional lung segmentation
  • Pulmonary artery segmentation
  • Pulmonary airway segmentation
  • CheXNet

• Mouse applications

  • Mouse brain extraction
  • Mouse brain parcellation
  • Mouse cortical thickness

• General applications

  • MRI super resolution
  • No reference image quality assesment using TID
  • Full reference image quality assessment

• Data augmentation

  • Noise
  • Histogram intensity warping
  • Simulate bias field
  • Random spatial transformations
  • Combined

---

b_0_ants_interfaces_basic.md

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

B-spline fitting

Python

>>> import ants
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> 
>>> x = np.linspace(-4, 4, num=100)
>>> y = np.exp(-np.multiply(x, x)) + np.random.uniform(-0.1, 0.1, len(x))
>>> u = np.linspace(0, 1.0, num=len(x))
>>> scattered_data = np.column_stack((x, y))
>>> parametric_data = np.expand_dims(u, axis=-1)
>>> spacing = 1/(len(x)-1)
>>> bspline_curve = ants.fit_bspline_object_to_scattered_data(scattered_data, 
...   parametric_data,
...   parametric_domain_origin=[0.0], parametric_domain_spacing=[spacing],
...   parametric_domain_size=[len(x)], is_parametric_dimension_closed=None,
...   number_of_fitting_levels=5, mesh_size=1)
>>> plt.plot(x, y, 'o', label='Noisy points')
[<matplotlib.lines.Line2D object at 0x31b572810>]
>>> plt.plot(bspline_curve[:,0], bspline_curve[:,1], label='B-spline curve')
[<matplotlib.lines.Line2D object at 0x31b53b200>]
>>> plt.grid(True)
>>> plt.axis('tight')
(-4.4, 4.4, -0.1478762733173202, 1.02880490107376)
>>> plt.legend(loc='upper left')
<matplotlib.legend.Legend object at 0x31b235340>
>>> plt.show()
2025-02-17 10:16:47.919 python3[17382:170684] +[IMKClient subclass]: chose IMKClient_Modern
2025-02-17 10:16:47.919 python3[17382:170684] +[IMKInputSession subclass]: chose IMKInputSession_Modern
>>>
>>> # 2-D closed curve example (hypocycloid)
>>> 
>>> number_of_sample_points = 100
>>> delta = 1/(number_of_sample_points)
>>> u = np.linspace(0, 1.0-delta, num=number_of_sample_points)
>>> x = 9 * np.cos(u * 2 * np.pi) + 3 * np.cos(3 * u * 2 * np.pi) + np.random.uniform(-1, 1, len(u))
>>> y = 9 * np.sin(u * 2 * np.pi) - 3 * np.sin(3 * u * 2 * np.pi) + np.random.uniform(-1, 1, len(u))
>>> scattered_data = np.column_stack((x, y))
>>> parametric_data = np.expand_dims(u, axis=-1)
>>> spacing = 1/(len(x)-1)
>>> bspline_curve = ants.fit_bspline_object_to_scattered_data(scattered_data, 
...   parametric_data,
...   parametric_domain_origin=[0.0], parametric_domain_spacing=[spacing],
...   parametric_domain_size=[len(x)], is_parametric_dimension_closed=[True],
...   number_of_fitting_levels=3, mesh_size=4)
>>> plt.plot(x, y, 'o', label='Noisy points')
[<matplotlib.lines.Line2D object at 0x32484d0d0>]
>>> plt.plot(bspline_curve[:,0], bspline_curve[:,1], label='B-spline curve')
[<matplotlib.lines.Line2D object at 0x3248b1550>]
>>> plt.grid(True)
>>> plt.axis('tight')
(-14.009181683424856, 14.060594273893225, -13.874089327839798, 12.798694412376495)
>>> plt.legend(loc='upper left')
<matplotlib.legend.Legend object at 0x32483fd70>
>>> plt.show()
>>>
>>> # 3-D B-spline surface of a trefoil knot
>>> 
>>> import ants
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> 
>>> uv = np.meshgrid(np.linspace(-2.0 * np.pi, 2.0 * np.pi, 100), 
...                  np.linspace(-1.0 * np.pi, 1.0 * np.pi, 100))
>>> uR = uv[0].flatten()
>>> vR = uv[1].flatten()
>>> 
>>> u = (uR + 2.0 * np.pi) / (4.0 * np.pi)
>>> v = (vR + 1.0 * np.pi) / (2.0 * np.pi)
>>> spacing = 1/99
>>> 
>>> x = np.cos(uR) * np.cos(vR) + 3.0 * np.cos(uR) * (1.5 + 0.5 * np.sin(1.5 * uR))
>>> y = np.sin(uR) * np.cos(vR) + 3.0 * np.sin(uR) * (1.5 + 0.5 * np.sin(1.5 * uR))
>>> z = np.sin(vR) + 2.0 * np.cos(1.5 * uR)
>>> 
>>> scattered_data_x = np.expand_dims(x, axis=-1)
>>> scattered_data_y = np.expand_dims(y, axis=-1)
>>> scattered_data_z = np.expand_dims(z, axis=-1)
>>> 
>>> parametric_data = np.column_stack((u, v))
>>> 
>>> bspline_mesh_x = ants.fit_bspline_object_to_scattered_data(
...   scattered_data_x, parametric_data,
...   parametric_domain_origin=[0.0, 0.0], parametric_domain_spacing=[0.1 * spacing, spacing],
...   parametric_domain_size=[1000, 100], is_parametric_dimension_closed=[True, True],
...   number_of_fitting_levels=4, mesh_size=[4, 7])
>>> bspline_mesh_y = ants.fit_bspline_object_to_scattered_data(
...   scattered_data_y, parametric_data,
...   parametric_domain_origin=[0.0, 0.0], parametric_domain_spacing=[0.1 * spacing, spacing],
...   parametric_domain_size=[1000, 100], is_parametric_dimension_closed=[True, True],
...   number_of_fitting_levels=4, mesh_size=[4, 7])
>>> bspline_mesh_z = ants.fit_bspline_object_to_scattered_data(
...   scattered_data_z, parametric_data,
...   parametric_domain_origin=[0.0, 0.0], parametric_domain_spacing=[0.1 * spacing, spacing],
...   parametric_domain_size=[1000, 100], is_parametric_dimension_closed=[True, True],
...   number_of_fitting_levels=4, mesh_size=[4, 7])
>>> 
>>> xx = bspline_mesh_x.numpy()
>>> yy = bspline_mesh_y.numpy()
>>> zz = bspline_mesh_z.numpy()
>>> 
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111, projection='3d')
>>> ax.plot_surface(xx, yy, zz)
<mpl_toolkits.mplot3d.art3d.Poly3DCollection object at 0x3248b26c0>
>>> plt.show()
>>> 
>>> # 2-D B-spline image approximation
>>> 
>>> number_of_random_points = 10000
>>> img = ants.image_read( ants.get_ants_data("r16"))
>>> img_array = img.numpy()
>>> row_indices = np.random.choice(range(2, img_array.shape[0]), number_of_random_points)
>>> col_indices = np.random.choice(range(2, img_array.shape[1]), number_of_random_points)
>>> scattered_data = np.zeros((number_of_random_points, 1))
>>> parametric_data = np.zeros((number_of_random_points, 2))
>>> 
>>> for i in range(number_of_random_points):
...     scattered_data[i,0] = img_array[row_indices[i], col_indices[i]]
...     parametric_data[i,0] = row_indices[i]
...     parametric_data[i,1] = col_indices[i]
... 
>>> bspline_img = ants.fit_bspline_object_to_scattered_data(
...     scattered_data, parametric_data,
...     parametric_domain_origin=[0.0, 0.0],
...     parametric_domain_spacing=[1.0, 1.0],
...     parametric_domain_size = img.shape,
...     number_of_fitting_levels=7, mesh_size=1)
>>> 
>>> f, ax = plt.subplots(1, 2)
>>> ax[0].imshow(np.rot90(img.numpy(), 1), cmap='gray')
<matplotlib.image.AxesImage object at 0x3672018e0>
>>> ax[0].set_title('Original image')
Text(0.5, 1.0, 'Original image')
>>> ax[0].set_xticklabels([])
[Text(-100.0, 0, ''), Text(0.0, 0, ''), Text(100.0, 0, ''), Text(200.0, 0, ''), Text(300.0, 0, '')]
>>> ax[0].set_yticklabels([])
[Text(0, -50.0, ''), Text(0, 0.0, ''), Text(0, 50.0, ''), Text(0, 100.0, ''), Text(0, 150.0, ''), Text(0, 200.0, ''), Text(0, 250.0, ''), Text(0, 300.0, '')]
>>> ax[1].imshow(np.rot90(bspline_img.numpy(), 1), cmap='gray')
<matplotlib.image.AxesImage object at 0x36725c320>
>>> ax[1].set_title('B-spline image')
Text(0.5, 1.0, 'B-spline image')
>>> ax[1].set_xticklabels([])
[Text(-100.0, 0, ''), Text(0.0, 0, ''), Text(100.0, 0, ''), Text(200.0, 0, ''), Text(300.0, 0, '')]
>>> ax[1].set_yticklabels([])
[Text(0, -50.0, ''), Text(0, 0.0, ''), Text(0, 50.0, ''), Text(0, 100.0, ''), Text(0, 150.0, ''), Text(0, 200.0, ''), Text(0, 250.0, ''), Text(0, 300.0, '')]
>>> plt.show()

R

> # 2-D curve example
> 
> x <- seq( from = -4, to = 4, by = 0.1 )
> y <- exp( -( x * x ) ) + runif( length( x ), min = -0.1, max = 0.1 )
> u <- seq( from = 0.0, to = 1.0, length.out = length( x ) )
> scatteredData <- cbind( x, y )
> parametricData <- as.matrix( u, ncol = 1 )
> numberOfSamplePoints <- 100
> spacing <- 1 / ( numberOfSamplePoints - 1 )
> bsplineCurve <- fitBsplineObjectToScatteredData( scatteredData, parametricData,
+    parametricDomainOrigin = c( 0.0 ), parametricDomainSpacing = c( spacing ),
+    parametricDomainSize = c( numberOfSamplePoints ), isParametricDimensionClosed = c( FALSE ),
+    numberOfFittingLevels = 5, meshSize = 1 )
> plot( x, y, "p", col = "red" )
> points( scatteredData[, 1], scatteredData[, 2], col = "green" )
> lines( bsplineCurve[, 1], bsplineCurve[, 2], col = "blue" )
> legend( "topleft", legend = c( "Noisy points", "B-spline curve" ), 
+        pch = c( 'o', '-' ), col = c( "green", "blue" ) )
> 
> # 2-D closed curve example (hypocycloid)
> 
> numberOfSamplePoints <- 100
> delta <- 1 / numberOfSamplePoints
> u <- seq(from = 0.0, to = 1.0 - delta, by = delta)
> x <- 9 * cos(u * 2 * pi) + 3 * cos(3 * u * 2 * pi) + runif( length( u ), min = -1, max = 1 )
> y <- 9 * sin(u * 2 * pi) - 3 * sin(3 * u * 2 * pi) + runif( length( u ), min = -1, max = 1 )
> scatteredData <- cbind( x, y )
> parametricData <- as.matrix( u, ncol = 1 )
> spacing <- 1 / ( numberOfSamplePoints - 1 )
> bsplineCurve <- fitBsplineObjectToScatteredData(scatteredData, parametricData,
+   parametricDomainOrigin = c( 0.0 ), parametricDomainSpacing = c( spacing ),
+   parametricDomainSize = c( numberOfSamplePoints ), isParametricDimensionClosed = c( TRUE ),
+   numberOfFittingLevels = 3, meshSize = 4 )
> plot( x, y, "p", col = "red" )
> points( scatteredData[, 1], scatteredData[, 2], col = "green" )
> lines( bsplineCurve[, 1], bsplineCurve[, 2], col = "blue" )
> legend( "topleft", legend = c("Noisy points", "B-spline-curve" ), 
+        pch = c( 'o', '-' ), col = c( "green", "blue" ) )
> 
> # 3-D B-spline surface of a trefoil knot
> 
> library( plot3D )
> 
> meshGrid <- function( u, v )
+ {
+   uv <- as.vector( t( outer( u, rep( 1, length( v ) ) ) ) )
+   vu <- as.vector( t( outer( rep( 1, length( u ) ), v ) ) )
+   return( list( uv, vu ) )
+ }
> 
> uv <- meshGrid( seq( from = -2.0 * pi, to = 2.0 * pi, by = 0.1 ), 
+                 seq( from = -1.0 * pi, to = 1.0 * pi, by = 0.1 ) )
> uR <- uv[[1]]
> vR <- uv[[2]]
> 
> u <- ( uR + 2.0 * pi ) / ( 4.0 * pi )
> v <- ( vR + 1.0 * pi ) / ( 2.0 * pi )
> 
> x <- cos( uR ) * cos( vR ) + 3.0 * cos( uR ) * ( 1.5 + 0.5 * sin( 1.5 * uR ) )
> y <- sin( uR ) * cos( vR ) + 3.0 * sin( uR ) * ( 1.5 + 0.5 * sin( 1.5 * uR ) )
> z <- sin( vR ) + 2.0 * cos( 1.5 * uR )
> 
> scatteredDataX <- as.matrix( x, ncol = 1 )
> scatteredDataY <- as.matrix( y, ncol = 1 )
> scatteredDataZ <- as.matrix( z, ncol = 1 )
> 
> parametricData <- cbind( u, v )
> 
> bsplineMeshX <- fitBsplineObjectToScatteredData(scatteredDataX, parametricData,
+   parametricDomainOrigin = c( 0.0, 0.0 ), parametricDomainSpacing = c( 0.001, 0.01 ),
+   parametricDomainSize = c( 1000, 100 ), isParametricDimensionClosed = c( TRUE, TRUE ),
+   numberOfFittingLevels = 4, meshSize = c( 4, 7 ) )
> bsplineMeshY <- fitBsplineObjectToScatteredData(scatteredDataY, parametricData,
+   parametricDomainOrigin = c( 0.0, 0.0 ), parametricDomainSpacing = c( 0.001, 0.01 ),
+   parametricDomainSize = c( 1000, 100 ), isParametricDimensionClosed = c( TRUE, TRUE ),
+   numberOfFittingLevels = 4, meshSize = c( 4, 7 ) )
> bsplineMeshZ <- fitBsplineObjectToScatteredData(scatteredDataZ, parametricData,
+   parametricDomainOrigin = c( 0.0, 0.0 ), parametricDomainSpacing = c( 0.001, 0.01 ),
+   parametricDomainSize = c( 1000, 100 ), isParametricDimensionClosed = c( TRUE, TRUE ),
+   numberOfFittingLevels = 4, meshSize = c( 4, 7 ) )
> 
> xx <- as.array( bsplineMeshX )
> yy <- as.array( bsplineMeshY )
> zz <- as.array( bsplineMeshZ )
> 
> surf3D( x = xx, y = yy, z = zz )
> 
> # 2-D B-spline image approximation
> 
> numberOfRandomPoints <- 10000
> img <- antsImageRead( getANTsRData( "r16" ) )
> imgArray <- as.array( img )
> rowIndices <- sample( 2:( dim( imgArray)[1] - 1 ), numberOfRandomPoints, replace = TRUE )
> colIndices <- sample( 2:( dim( imgArray)[2] - 1 ), numberOfRandomPoints, replace = TRUE )
> 
> scatteredData <- as.matrix( array( data = 0, dim = c( numberOfRandomPoints, 1 ) ) )
> parametricData <- as.matrix( array( data = 0, dim = c( numberOfRandomPoints, 2 ) ) )
> for( i in seq_len( numberOfRandomPoints ) )
+   {
+   scatteredData[i, 1] <- imgArray[rowIndices[i], colIndices[i]]
+   parametricData[i, 1] <- rowIndices[i]
+   parametricData[i, 2] <- colIndices[i]
+   }
> 
> bsplineImage <- fitBsplineObjectToScatteredData( scatteredData, parametricData,
+   parametricDomainOrigin = c( 0.0, 0.0 ), parametricDomainSpacing = c( 1.0, 1.0 ),
+   parametricDomainSize = dim( img ),
+   numberOfFittingLevels = 7, meshSize = c( 1, 1 ) )
> 
> img <- iMath( img, operation = "Normalize" )
> bsplineImage <- iMath( bsplineImage, operation = "Normalize" )
> 
> nf <- layout( matrix( c( 1, 2 ), ncol = 2 ), 1, 1 )
> image( as.array( img ), col = gray.colors( 100, start = 0.0, end = 1.0 ), axes = FALSE )
> title( main = "Original image", font.main = 4 )
> image( as.array( bsplineImage ), col = gray.colors( 100, start = 0.0, end = 1.0 ), axes = FALSE )
> title( main = "B-spline image", font.main = 4 )

b_1_ants_interfaces_image_registration.md

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 ) 

Label image registration

Python

>>> import ants
>>> import antspynet
>>> 
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r')) 
>>> dkt = antspynet.desikan_killiany_tourville_labeling(t1, version=1)
>>> hoa = antspynet.harvard_oxford_atlas_labeling(t1)['segmentation_image']
>>> 
>>> t1_template = ants.image_read(antspynet.get_antsxnet_data('adni')) 
>>> dkt_template = antspynet.desikan_killiany_tourville_labeling(t1_template, version=1)
>>> hoa_template = antspynet.harvard_oxford_atlas_labeling(t1_template)['segmentation_image']
>>> 
>>> reg = ants.label_image_registration([dkt_template, hoa_template],
                                        [dkt, hoa],  
                                        fixed_intensity_images=t1_template,
                                        moving_intensity_images=t1,
                                        initial_transforms='affine',
                                        type_of_deformable_transform='antsRegistrationSyNQuick[bo]',
                                        label_image_weighting=[2.0, 1.0],
                                        verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
> 
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) ) 
> dkt <- desikanKillianyTourvilleLabeling( t1, version = 1 )
> hoa <- harvardOxfordAtlasLabeling( t1 )$segmentationImage 
> 
> t1Template <- antsImageRead( getANTsXNetData( 'adni' ) ) 
> dktTemplate <- desikanKillianyTourvilleLabeling( t1Template, version = 1 )
> hoaTemplate <- harvardOxfordAtlasLabeling( t1Template )$segmentationImage 
> 
> reg <- labelImageRegistration(list( dktTemplate, hoaTemplate ),
                                list( dkt, hoa ),  
                                fixedIntensityImages = t1Template,
                                movingIntensityImages = t1,
                                initialTransforms = 'affine',
                                typeOfDeformableTransform = 'antsRegistrationSyNQuick[bo]',
                                labelImageWeighting = c( 2.0, 1.0 ),
                                verbose = TRUE )

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

b_2_ants_interfaces_landmark_registration.md

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 )

b_3_ants_velocity_flows.md

Velocity flows across point sets

Two more complicated examples:

• Developmental mouse brain mapping example
• 4D lung CT

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

b_4_ants_interfaces_misc_utilities.md

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

b_5_ants_interfaces_deformation_based_morphometry.md

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 )

c_0_antsxnet_basic.md

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

c_1_antsxnet_applications_1_brain.md

ANTsXNet

Brain applications

Brain extraction

• T1 (multiple "flavors")
• T1 (three-tissue)
• T1 (hemispheres)
• T1 (lobes)
• 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)
>>> 
>>> # Three-tissue
>>> bext = antspynet.brain_extraction(t1, modality="t1threetissue", verbose=True)
>>> seg = bext['segmentation_image']
>>> ants.plot(t1, overlay=seg, overlay_alpha=0.5)
>>>
>>> # Hemispheres
>>> bext = antspynet.brain_extraction(t1, modality="t1hemi", verbose=True)
>>> seg = bext['segmentation_image']
>>> ants.plot(t1, overlay=seg, overlay_alpha=0.5)
>>>
>>> # Lobes
>>> bext = antspynet.brain_extraction(t1, modality="t1lobes", verbose=True)
>>> seg = bext['segmentation_image']
>>> 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 )
>
> # Three-tissue
> bext <- brainExtraction( t1, modality = "t1threetissue", verbose = TRUE )
> seg <- bext$segmentationImage
> plot( t1, seg, alpha = 0.5 )
>
> # Hemispheres
> bext <- brainExtraction( t1, modality = "t1hemi", verbose = TRUE )
> seg <- bext$segmentationImage
> plot( t1, seg, alpha = 0.5 )
>
> # Lobes
> bext <- brainExtraction( t1, modality = "t1lobes", verbose = TRUE )
> seg <- bext$segmentationImage
> plot( t1, seg, alpha = 0.5 )

Deep Atropos

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 )

Cortical thickness

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(...)

Desikan-Killiany-Tourville parcellation

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('mprage_hippmapp3r'))
>>>
>>> # original version
>>> dkt = antspynet.desikan_killiany_tourville_labeling(t1, do_lobar_parcellation=True, version=0, verbose=True)
>>> # updated version
>>> dkt = antspynet.desikan_killiany_tourville_labeling(t1, do_lobar_parcellation=True, version=1, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
>
> # original version
> dkt <- desikanKillianyTourvilleLabeling( t1, doLobarParcellation = TRUE, version = 0, verbose = TRUE )
> # updated version
> dkt <- desikanKillianyTourvilleLabeling( t1, doLobarParcellation = TRUE, version = 1, verbose = TRUE )

Harvard-Oxford Atlas subcortical parcellation

Python

>>> import ants
>>> import antspynet
>>>
>>> t1 = ants.image_read(antspynet.get_antsxnet_data('adni'))
>>> hoa = antspynet.harvard_oxford_atlas_labeling(t1, verbose=True)

R

> library( ANTsR )
> library( ANTsRNet )
> 
> t1 <- antsImageRead( getANTsXNetData( "adni" ) )
> hoa <- harvardOxfordAtlasLabeling( t1, verbose = TRUE )

DeepFLASH (medial temporal lobe parcellation)

• 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 )

Hippmapp3r (hippocampal segmentation)

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 )

Brain Age

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 )

Claustrum segmentation

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 )

Hypothalamus segmentation

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 )

White matter hyperintensities segmentation (SYSU)

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 )

White matter hyperintensities segmentation (Hypermapp3r)

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 )

White matter hyperintensities segmentation (SHIVA)

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 )

White matter hyperintensities segmentation (ANTsXNet)

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 )

Perivascular spaces segmentation (SHIVA)

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" )	

Brain tumor segmentation

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" )

MRA-TOF vessel segmentation

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 )

Lesion segmentation (WIP)

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" )

Whole head inpainting

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" )

c_1_antsxnet_applications_2_lung.md

Lung applications

Lung extraction

• 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 )

Functional lung segmentation

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 )

Pulmonary artery segmentation

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 )

Pulmonary airway segmentation

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 )

CheXNet

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

---

c_1_antsxnet_applications_3_mouse.md

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 )

c_1_antsxnet_applications_4_general.md

General applications

MRI super resolution

Python

>>> import ants
>>> import antspynet
>>> 
>>> 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, expansion_factor=[1,2,2], verbose=True)
>>> ants.plot(t1_sr)

R

> library( ANTsR )
> library( ANTsRNet )
>
> t1 <- antsImageRead( getANTsXNetData( 'mprage_hippmapp3r' ) )
> t1Lr <- resampleImage( t1, c( 4, 4, 4 ), useVoxels = FALSE )
> plot( t1Lr )
> t1Sr <- mriSuperResolution( t1Lr, expansionFactor = c( 1, 2, 2 ), 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 ) )

c_2_antsxnet_data_augmentation.md

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" ) )
>   }