<jupyter_start><jupyter_text>SimpleITK Spatial Transformations**Summary:**1. Points are represented by vector-like data types: Tuple, Numpy array, List.2. Matrices are represented by vector-like data types in row major order.3. Default transformation initialization as the identity transform.4. Angles specified in radians, distances specified in unknown but consistent units (nm,mm,m,km...).5. All global transformations **except translation** are of the form:$$T(\mathbf{x}) = A(\mathbf{x}-\mathbf{c}) + \mathbf{t} + \mathbf{c}$$   Nomenclature (when printing your transformation):   * Matrix: the matrix $A$   * Center: the point $\mathbf{c}$   * Translation: the vector $\mathbf{t}$   * Offset: $\mathbf{t} + \mathbf{c} - A\mathbf{c}$6. Bounded transformations, BSplineTransform and DisplacementFieldTransform, behave as the identity transform outside the defined bounds.7. DisplacementFieldTransform:   * Initializing the DisplacementFieldTransform using an image requires that the image's pixel type be sitk.sitkVectorFloat64.   * Initializing the DisplacementFieldTransform using an image will "clear out" your image (your alias to the image will point to an empty, zero sized, image).8. Composite transformations are applied in stack order (first added, last applied).  Transformation TypesSimpleITK supports the following transformation types.TranslationTransform2D or 3D, translation  VersorTransform3D, rotation represented by a versor  VersorRigid3DTransform3D, rigid transformation with rotation represented by a versor  Euler2DTransform2D, rigid transformation with rotation represented by a Euler angle  Euler3DTransform3D, rigid transformation with rotation represented by Euler angles  Similarity2DTransform2D, composition of isotropic scaling and rigid transformation with rotation represented by a Euler angle  Similarity3DTransform3D, composition of isotropic scaling and rigid transformation with rotation represented by a versor  ScaleTransform2D or 3D, anisotropic scaling  ScaleVersor3DTransform3D, rigid transformation and anisotropic scale is added to the rotation matrix part (not composed as one would expect)  ScaleSkewVersor3DTransform3D, rigid transformation with anisotropic scale and skew matrices added to the rotation matrix part (not composed as one would expect)  AffineTransform2D or 3D, affine transformation.  BSplineTransform2D or 3D, deformable transformation represented by a sparse regular grid of control points.   DisplacementFieldTransform2D or 3D, deformable transformation represented as a dense regular grid of vectors.  Transform  A generic transformation. Can represent any of the SimpleITK transformations, and a composite transformation (stack of transformations concatenated via composition, last added, first applied).<jupyter_code>library(SimpleITK)

source("utilities.R")
OUTPUT_DIR <- "Output"<jupyter_output><empty_output><jupyter_text>We will introduce the transformation types, starting with translation and illustrating how to move from a lower to higher parameter space (e.g. translation to rigid).  We start with the global transformations. All of them except translation are of the form:$$T(\mathbf{x}) = A(\mathbf{x}-\mathbf{c}) + \mathbf{t} + \mathbf{c}$$In ITK speak (when printing your transformation):Matrix: the matrix $A$Center: the point $\mathbf{c}$Translation: the vector $\mathbf{t}$Offset: $\mathbf{t} + \mathbf{c} - A\mathbf{c}$  TranslationTransformCreate a translation and then transform a point and use the inverse transformation to get the original back.<jupyter_code>dimension <- 2        
offset <- rep(1,dimension)*2 # Use a "trick" to create a vector of 2s based on dimension. 
translation <- TranslationTransform(dimension, offset)
print(translation)
point <- if(dimension==2) c(10, 11) else c(10, 11, 12) # Set point to match dimension
transformed_point <- translation$TransformPoint(point)
translation_inverse <- translation$GetInverse()
cat(paste0("original point: ", point2str(point), "\n",
          "transformed point: ", point2str(transformed_point), "\n",
          "back to original: ", point2str(translation_inverse$TransformPoint(transformed_point))))<jupyter_output><empty_output><jupyter_text>Euler2DTransformRigidly transform a 2D point using a Euler angle parameter specification.Notice that the dimensionality of the Euler angle based rigid transformation is associated with the class, unlike the translation which is set at construction.<jupyter_code>point <- c(10, 11)
rotation2D <- Euler2DTransform()
rotation2D$SetTranslation(c(7.2, 8.4))
rotation2D$SetAngle(pi/2.0)
cat(paste0("original point: ", point2str(point), "\n",
      "transformed point: ", point2str(rotation2D$TransformPoint(point)),"\n"))<jupyter_output><empty_output><jupyter_text>VersorTransform (rotation in 3D)Rotation using a versor, vector part of unit quaternion, parameterization. Quaternion defined by rotation of $\theta$ radians around axis $n$, is $q = [n*\sin(\frac{\theta}{2}), \cos(\frac{\theta}{2})]$.<jupyter_code># Use a versor:
rotation1 <- VersorTransform(c(0,0,1,0))

# Use axis-angle:
rotation2 <- VersorTransform(c(0,0,1), pi)

# Use a matrix:
rotation3 <- VersorTransform()
rotation3$SetMatrix(c(-1, 0, 0, 0, -1, 0, 0, 0, 1))

point <- c(10, 100, 1000)

p1 <- rotation1$TransformPoint(point)
p2 <- rotation2$TransformPoint(point)
p3 <- rotation3$TransformPoint(point)

cat(paste0("Points after transformation:\np1=", point2str(p1,15), 
      "\np2=", point2str(p2,15),"\np3=", point2str(p3,15)))<jupyter_output><empty_output><jupyter_text>Translation to Rigid [3D]We only need to copy the translational component.<jupyter_code>dimension <- 3        
trans <- c(1,2,3) 
translation <- TranslationTransform(dimension, trans)

# Copy the translational component.
rigid_euler <- Euler3DTransform()
rigid_euler$SetTranslation(translation$GetOffset())

# Apply the transformations to the same set of random points and compare the results.
print_transformation_differences(translation, rigid_euler)<jupyter_output><empty_output><jupyter_text>Rotation to Rigid [3D]Copy the matrix or versor and center of rotation.<jupyter_code>rotationCenter <- c(10, 10, 10)
rotation <- VersorTransform(c(0,0,1,0), rotationCenter)

rigid_versor <- VersorRigid3DTransform()
rigid_versor$SetRotation(rotation$GetVersor())
#rigid_versor$SetCenter(rotation$GetCenter()) # Intentional error, not copying center of rotation

# Apply the transformations to the same set of random points and compare the results.
print_transformation_differences(rotation, rigid_versor)<jupyter_output><empty_output><jupyter_text>In the cell above, when we don't copy the center of rotation we have a constant error vector, $\mathbf{c}$ - A$\mathbf{c}$.  Similarity [2D]When the center of the similarity transformation is not at the origin the effect of the transformation is not what most of us expect. This is readily visible if we limit the transformation to scaling: $T(\mathbf{x}) = s\mathbf{x}-s\mathbf{c} + \mathbf{c}$. Changing the transformation's center results in scale + translation.<jupyter_code># 2D square centered on (0,0)
points <- matrix(data=c(-1.0,-1.0, -1.0,1.0, 1.0,1.0, 1.0,-1.0), ncol=4, nrow=2) 
# Scale by 2 (center default is [0,0])
similarity <- Similarity2DTransform();
similarity$SetScale(2)

scaled_points <- apply(points, MARGIN=2, similarity$TransformPoint) 

#Uncomment the following lines to change the transformations center and see what happens:
#similarity$SetCenter(c(0,2))
#scaled_points <- apply(points, 2, similarity$TransformPoint) 

plot(points[1,],points[2,], xlim=c(-10,10), ylim=c(-10,10), pch=19, col="blue", xlab="", ylab="", las=1)
points(scaled_points[1,], scaled_points[2,], col="red", pch=17)
legend('top', col= c("red", "blue"), pch=c(17,19), legend = c("transformed points", "original points"))<jupyter_output><empty_output><jupyter_text>Rigid to Similarity [3D]Copy the translation, center, and matrix or versor.<jupyter_code>rotation_center <- c(100, 100, 100)
theta_x <- 0.0
theta_y <- 0.0
theta_z <- pi/2.0
translation <- c(1,2,3)

rigid_euler <- Euler3DTransform(rotation_center, theta_x, theta_y, theta_z, translation)

similarity <- Similarity3DTransform()
similarity$SetMatrix(rigid_euler$GetMatrix())
similarity$SetTranslation(rigid_euler$GetTranslation())
similarity$SetCenter(rigid_euler$GetCenter())

# Apply the transformations to the same set of random points and compare the results.
print_transformation_differences(rigid_euler, similarity)<jupyter_output><empty_output><jupyter_text>Similarity to Affine [3D]Copy the translation, center and matrix.<jupyter_code>rotation_center <- c(100, 100, 100)
axis <- c(0,0,1)
angle <- pi/2.0
translation <- c(1,2,3)
scale_factor <- 2.0
similarity <- Similarity3DTransform(scale_factor, axis, angle, translation, rotation_center)

affine <- AffineTransform(3)
affine$SetMatrix(similarity$GetMatrix())
affine$SetTranslation(similarity$GetTranslation())
affine$SetCenter(similarity$GetCenter())

# Apply the transformations to the same set of random points and compare the results.
print_transformation_differences(similarity, affine)<jupyter_output><empty_output><jupyter_text>Scale TransformJust as the case was for the similarity transformation above, when the transformations center is not at the origin, instead of a pure anisotropic scaling we also have translation ($T(\mathbf{x}) = \mathbf{s}^T\mathbf{x}-\mathbf{s}^T\mathbf{c} + \mathbf{c}$).<jupyter_code># 2D square centered on (0,0).
points <- matrix(data=c(-1.0,-1.0, -1.0,1.0, 1.0,1.0, 1.0,-1.0), ncol=4, nrow=2) 

# Scale by half in x and 2 in y.
scale <- ScaleTransform(2, c(0.5,2));

scaled_points <- apply(points, 2, scale$TransformPoint) 

#Uncomment the following lines to change the transformations center and see what happens:
#scale$SetCenter(c(0,2))
#scaled_points <- apply(points, 2, scale$TransformPoint) 

plot(points[1,],points[2,], xlim=c(-10,10), ylim=c(-10,10), pch=19, col="blue", xlab="", ylab="", las=1)
points(scaled_points[1,], scaled_points[2,], col="red", pch=17)
legend('top', col= c("red", "blue"), pch=c(17,19), legend = c("transformed points", "original points"))<jupyter_output><empty_output><jupyter_text>Unintentional Misnomers (originally from ITK)Two transformation types whose names may mislead you are ScaleVersor and ScaleSkewVersor. Basing your choices on expectations without reading the documentation will surprise you.ScaleVersor -  based on name expected a composition of transformations, in practice it is:$$T(x) = (R+S)(\mathbf{x}-\mathbf{c}) + \mathbf{t} + \mathbf{c},\;\; \textrm{where } S= \left[\begin{array}{ccc} s_0-1 & 0 & 0 \\ 0 & s_1-1 & 0 \\ 0 & 0 & s_2-1 \end{array}\right]$$ ScaleSkewVersor - based on name expected a composition of transformations, in practice it is:$$T(x) = (R+S+K)(\mathbf{x}-\mathbf{c}) + \mathbf{t} + \mathbf{c},\;\; \textrm{where } S = \left[\begin{array}{ccc} s_0-1 & 0 & 0 \\ 0 & s_1-1 & 0 \\ 0 & 0 & s_2-1 \end{array}\right]\;\; \textrm{and } K = \left[\begin{array}{ccc} 0 & k_0 & k_1 \\ k_2 & 0 & k_3 \\ k_4 & k_5 & 0 \end{array}\right]$$ Note that ScaleSkewVersor is is an over-parametrized version of the affine transform, 15 parameters (scale, skew, versor, translation) vs. 12 parameters (matrix, translation).  Bounded TransformationsSimpleITK supports two types of bounded non-rigid transformations, BSplineTransform (sparse representation) and 	DisplacementFieldTransform (dense representation).Transforming a point that is outside the bounds will return the original point - identity transform.  BSplineUsing a sparse set of control points to control a free form deformation. Using the cell below it is clear that the BSplineTransform allows for folding and tearing.<jupyter_code># Create the transformation (when working with images it is easier to use the BSplineTransformInitializer function
# or its object oriented counterpart BSplineTransformInitializerFilter).
dimension <- 2
spline_order <- 3
direction_matrix_row_major <- c(1.0,0.0,0.0,1.0) # identity, mesh is axis aligned
origin <- c(-1.0,-1.0)  
domain_physical_dimensions <- c(2,2)

bspline <- BSplineTransform(dimension, spline_order)
bspline$SetTransformDomainOrigin(origin)
bspline$SetTransformDomainDirection(direction_matrix_row_major)
bspline$SetTransformDomainPhysicalDimensions(domain_physical_dimensions)
bspline$SetTransformDomainMeshSize(c(4,3))

# Random displacement of the control points.
originalControlPointDisplacements <- runif(length(bspline$GetParameters()))
bspline$SetParameters(originalControlPointDisplacements)

# Apply the BSpline transformation to a grid of points 
# starting the point set exactly at the origin of the BSpline mesh is problematic as
# these points are considered outside the transformation's domain,
# remove epsilon below and see what happens.
numSamplesX <- 10
numSamplesY <- 20

eps <- .Machine$double.eps

coordsX <- seq(origin[1] + eps,
               origin[1] + domain_physical_dimensions[1],
               (domain_physical_dimensions[1]-eps)/(numSamplesX-1))
coordsY <- seq(origin[2] + eps,
               origin[2] + domain_physical_dimensions[2],
               (domain_physical_dimensions[2]-eps)/(numSamplesY-1))
# next two lines equivalent to Python's/MATLAB's meshgrid 
XX <- outer(coordsY*0, coordsX, "+")
YY <- outer(coordsY, coordsX*0, "+")  

display_displacement_scaling_effect(0.0, XX, YY, bspline, originalControlPointDisplacements)

# Visualize the effect of scaling the control point displacements 
# on our set of points (we recommend keeping the scaling in the 
# range [-1.5,1.5] due to display bounds). 
display_displacement_scaling_effect(s=0.5, XX, YY, bspline, originalControlPointDisplacements)<jupyter_output><empty_output><jupyter_text>DisplacementFieldA dense set of vectors representing the displacement inside the given domain. The most generic representation of a transformation.<jupyter_code># Create the displacement field. 
    
# When working with images the safer thing to do is use the image based constructor,
# sitk.DisplacementFieldTransform(my_image), all the fixed parameters will be set correctly and the displacement
# field is initialized using the vectors stored in the image. SimpleITK requires that the image's pixel type be 
# sitk.sitkVectorFloat64.

displacement <- DisplacementFieldTransform(2)
field_size <- c(10,20)
field_origin <- c(-1.0,-1.0)  
field_spacing <- c(2.0/9.0,2.0/19.0)   
field_direction <- c(1,0,0,1) # direction cosine matrix (row major order)     

# Concatenate all the information into a single list
displacement$SetFixedParameters(c(field_size, field_origin, field_spacing, field_direction))
# Set the interpolator, either sitkLinear which is default or nearest neighbor
displacement$SetInterpolator("sitkNearestNeighbor")

originalDisplacements <- runif(length(displacement$GetParameters()))
displacement$SetParameters(originalDisplacements)

coordsX <- seq(field_origin[1],
               field_origin[1]+(field_size[1]-1)*field_spacing[1],
               field_spacing[1])
coordsY <- seq(field_origin[2],
               field_origin[2]+(field_size[2]-1)*field_spacing[2],
               field_spacing[2])

# next two lines equivalent to Python's/MATLAB's meshgrid 
XX <- outer(coordsY*0, coordsX, "+")
YY <- outer(coordsY, coordsX*0, "+")  

display_displacement_scaling_effect(0.0, XX, YY, displacement, originalDisplacements)

# Visualize the effect of scaling the control point displacements 
# on our set of points (we recommend keeping the scaling in the 
# range [-1.5,1.5] due to display bounds). 
display_displacement_scaling_effect(0.5, XX, YY, displacement, originalDisplacements)<jupyter_output><empty_output><jupyter_text>Composite transform (Transform)The generic SimpleITK transform class. This class can represent both a single transformation (global, local), or a composite transformation (multiple transformations applied one after the other). This is the output typed returned by the SimpleITK registration framework. The choice of whether to use a composite transformation or compose transformations on your own has subtle differences in the registration framework.Composite transforms enable a combination of a global transformation with multiple local/bounded transformations. This is useful if we want to apply deformations only in regions that deform while other regions are only effected by the global transformation.The following code illustrates this, where the whole region is translated and subregions have different deformations.<jupyter_code># Global transformation.
translation <- TranslationTransform(2, c(1.0,0.0))

# Displacement in region 1.
displacement1 <- DisplacementFieldTransform(2)
field_size <- c(10,20)
field_origin <- c(-1.0,-1.0)  
field_spacing <- c(2.0/9.0,2.0/19.0)   
field_direction <- c(1,0,0,1) # direction cosine matrix (row major order)     

# Concatenate all the information into  a single list.
displacement1$SetFixedParameters(c(field_size, field_origin, field_spacing, field_direction))
displacement1$SetParameters(rep(1.0, length(displacement1$GetParameters())))

# Displacement in region 2.
displacement2 <- DisplacementFieldTransform(2)
field_size <- c(10,20)
field_origin <- c(1.0,-3)  
field_spacing <- c(2.0/9.0,2.0/19.0)   
field_direction <- c(1,0,0,1) #direction cosine matrix (row major order)     

# Concatenate all the information into a single list.
displacement2$SetFixedParameters(c(field_size, field_origin, field_spacing, field_direction))
displacement2$SetParameters(rep(-1.0, length(displacement2$GetParameters())))

# Composite transform which applies the global and local transformations.
composite <- Transform(translation)
composite$AddTransform(displacement1)
composite$AddTransform(displacement2)

# Apply the composite transformation to points in ([-1,-3],[3,1]) and 
# display the deformation using a quiver plot.
        
# Generate points.
numSamplesX <- 10
numSamplesY <- 10
coordsX <- seq(-1.0, 3.0, 4.0/(numSamplesX-1))
coordsY <- seq(-3.0, 1.0, 4.0/(numSamplesY-1))
# next two lines equivalent to Python's/MATLAB's meshgrid 
original_x_mat <- outer(coordsY*0, coordsX, "+")
original_y_mat <- outer(coordsY, coordsX*0, "+")  

# Transform points and plot.
original_points <- mapply(function(x,y) c(x,y), original_x_mat, original_y_mat)
transformed_points <- mapply(function(x,y) composite$TransformPoint(c(x,y)), original_x_mat, original_y_mat)
plot(0,0,xlim=c(-1.0,3.0), ylim=c(-3.0,1.0), las=1)
arrows(original_points[1,], original_points[2,], transformed_points[1,], transformed_points[2,])<jupyter_output><empty_output><jupyter_text>Writing and ReadingThe SimpleITK.ReadTransform() returns a SimpleITK.Transform . The content of the file can be any of the SimpleITK transformations or a composite (set of transformations).<jupyter_code># Create a 2D rigid transformation, write it to disk and read it back.
basic_transform <- Euler2DTransform()
basic_transform$SetTranslation(c(1,2))
basic_transform$SetAngle(pi/2.0)

full_file_name <- file.path(OUTPUT_DIR, "euler2D.tfm")

WriteTransform(basic_transform, full_file_name)

# The ReadTransform function returns a SimpleITK Transform no matter the type of the transform 
# found in the file (global, bounded, composite).
read_result <- ReadTransform(full_file_name)
cat(paste("Original type: ",basic_transform$GetName(),"\nType after reading: ", read_result$GetName(),"\n"))
print_transformation_differences(basic_transform, read_result)


# Create a composite transform then write and read.
displacement <- DisplacementFieldTransform(2)
field_size <- c(10,20)
field_origin <- c(-10.0,-100.0)  
field_spacing <- c(20.0/(field_size[1]-1),200.0/(field_size[2]-1)) 
field_direction <- c(1,0,0,1) #direction cosine matrix (row major order)

# Concatenate all the information into a single list.
displacement$SetFixedParameters(c(field_size, field_origin, field_spacing, field_direction))
displacement$SetParameters(runif(length(displacement$GetParameters())))

composite_transform <- Transform(basic_transform)
composite_transform$AddTransform(displacement)

full_file_name <- file.path(OUTPUT_DIR, "composite.tfm")

WriteTransform(composite_transform, full_file_name)
read_result <- ReadTransform(full_file_name)
cat("\n")
print_transformation_differences(composite_transform, read_result)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10896170548756415, 'fraction_numerical': 0.027339123568890643, 'mean_word_length': 6.036458333333333, 'pattern_counts': {'":': 0, '<': 157, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '89928', 'license': 'Apache-2.0', 'n_tokens_mistral': 6519, 'n_tokens_neox': 6004, 'n_words': 1862, 'path': 'r/01_spatial_transformations.ipynb', 'repo_name': 'zivy/ISBI2018_TUTORIAL', 'size': '20264'}