Hi everyone. I’m new to SimpleITK. It’s been great so far and I’m hoping to do more with it.

I’m trying to perform a dense deformable registration between two 3D volumes, initialized with an affine registration performed prior (see `affine_matrix`

below). However, the moved image I’m getting contain artefacts at the edges of the image in all three dimensions (first and last planes, as well as the first and last pixels in x and y). I noticed that that if the shrink factors I set is not a factor of any of the size of the three dimensions, the remainder of the voxels would be set to a transformation of identity. For example, if there are 100 planes (i.e. z=100) and I set a shrink factor of 9, then one plane (the first plane) will be exactly the first plane of the moving image. This is true for the other dimensions as well.

Furthermore, if I set varianceForUpdateField to a non-zero number (`field_smoothing`

kwarg in `deformableAlign()`

below), the artefact returns (i.e. the first and last planes is identity or somewhere between identity and a somewhat accurate displacement vector in the other planes) even when the shrinkFactor is set to 1.

Here’s how I set up the deformable registration:

```
def deformableAlign(
fixed, moving,
fixed_vox=None,
moving_vox=None,
affine_matrix=None,
ncc_radius=8,
gradient_smoothing=1.0,
field_smoothing=0.5,
shrink_factors=[2,],
smooth_sigmas=[2,],
learning_rate=1.0,
number_of_iterations=250,
):
fixed, moving = numpyToSITK(fixed, moving, fixed_vox, moving_vox) # convert to sitk images, set spacing
affine = matrixToAffineTransform(affine_matrix)
irm = getDeformableRegistrationModel(
fixed_vox,
learning_rate,
number_of_iterations,
shrink_factors,
smooth_sigmas,
ncc_radius,
)
tdff = sitk.TransformToDisplacementFieldFilter()
tdff.SetReferenceImage(fixed)
df = tdff.Execute(affine)
dft = sitk.DisplacementFieldTransform(df) ## consume an image to construct a displacement field transform
dft.SetSmoothingGaussianOnUpdate(
varianceForUpdateField=gradient_smoothing,
varianceForTotalField=field_smoothing,
)
irm.SetInitialTransform(dft, inPlace=True)
deformation = irm.Execute(sitk.Cast(fixed, sitk.sitkFloat32),
sitk.Cast(moving, sitk.sitkFloat32),
)
tdff.SetOutputPixelType(sitk.sitkVectorFloat32)
return sitk.GetArrayFromImage(tdff.Execute(deformation))
def getDeformableRegistrationModel(
fixed_vox,
learning_rate,
iterations,
shrink_factors,
smooth_sigmas,
ncc_radius,
):
# set up registration object
ncores = int(os.environ["LSB_DJOB_NUMPROC"])
sitk.ProcessObject.SetGlobalDefaultNumberOfThreads(2*ncores)
irm = sitk.ImageRegistrationMethod()
irm.SetNumberOfThreads(2*ncores)
irm.SetInterpolator(sitk.sitkLinear)
# metric
irm.SetMetricAsANTSNeighborhoodCorrelation(ncc_radius) ## might take a list
irm.MetricUseFixedImageGradientFilterOff()
# optimizer
max_step = np.min(fixed_vox)
irm.SetOptimizerAsGradientDescent(
numberOfIterations=iterations,
learningRate=learning_rate,
maximumStepSizeInPhysicalUnits=max_step,
)
irm.SetOptimizerScalesFromPhysicalShift()
# pyramid
irm.SetShrinkFactorsPerLevel(shrinkFactors=shrink_factors)
irm.SetSmoothingSigmasPerLevel(smoothingSigmas=smooth_sigmas)
irm.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
# callback
def callback(irm):
level = irm.GetCurrentLevel()
iteration = irm.GetOptimizerIteration()
metric = irm.GetMetricValue()
print("LEVEL: ", level, " ITERATION: ", iteration, " METRIC: ", metric)
irm.AddCommand(sitk.sitkIterationEvent, lambda: callback(irm))
return irm
```

Any help is appreciated! Thanks!