## Landmark registration

### Python 

__Two simple examples__

```python
>>> 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__

```python
>>> 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__

```r
> 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__

```r
> 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 )
```