Why is my image registration code taking so long to run?

Question

I have written the following code for image registration for head CT images. Although it worked a few times previously, now it is taking very long time to finish running. In fact, it goes on for over 30 minutes without me interrupting the kernel. Even then, it will take a long time to respond to the interruption. I am running this code on Spyder.

Here is my code:

import SimpleITK as sitk
import numpy as np
import matplotlib.pyplot as plt


img_fix = sitk.ReadImage("template file for head ct") #template
img_fix_array = sitk.GetArrayFromImage(img_fix)
img_fix_size = img_fix.GetSize()
img_fix_spacing = img_fix.GetSpacing()

img = sitk.ReadImage("head ct scan image")
img_array = sitk.GetArrayFromImage(img)

desired_spacing = (0.4609375, 0.4609375, 0.5)
new_size = img.GetSize()

new_im_fix = sitk.Image(new_size, img_fix.GetPixelID())
new_im_fix.SetSpacing(desired_spacing)



# Resample the content of the original template image onto the new template image
resampler = sitk.ResampleImageFilter()
resampler.SetSize(new_size)
resampler.SetOutputSpacing(desired_spacing)
resampler.SetOutputOrigin(img_fix.GetOrigin())
resampler.SetOutputDirection(img_fix.GetDirection())
resampler.SetTransform(sitk.Transform())
new_im_fix = resampler.Execute(img_fix)


new_array = sitk.GetArrayFromImage(new_im_fix)

def _rescale(image, brain_only=False):  #what is brain_only ?
  """Rescale to match SCCT template intensities"""
  vol = sitk.GetArrayFromImage(img)
  vol_new = np.zeros(vol.shape)
  bw_air = vol <= -100
  vol_new[bw_air] = vol[bw_air] + 1000  # air: 0 to 900 (-Inf to -100)
  bw_brain = np.logical_and(vol > -100, vol <= 100)
  vol_new[bw_brain] = ((vol[bw_brain] + 99) / 199 * 2199) + 901  # brain: 901 to 3100
  bw_bone = vol > 100
  vol_new[bw_bone] = vol[bw_bone] + 3000  # bone: 3101 to Inf (101 to Inf)
  if brain_only:  # zero out non-brain tissue
     vol_new[bw_air] = 0
     vol_new[bw_bone] = 0
  img_new = sitk.GetImageFromArray(vol_new.astype(vol.dtype))
  img_new.SetOrigin(img.GetOrigin())  # copy over physical space metadata
  img_new.SetDirection(img.GetDirection())
  img_new.SetSpacing(img.GetSpacing())
  return img_new

 im_mov = _rescale(img, brain_only=True)

R = sitk.ImageRegistrationMethod()
# R.SetMetricFixedMask(bw_fix)
 # R.SetMetricMovingMask(bw_mov)
R.SetMetricAsANTSNeighborhoodCorrelation(1) #changed to ANTS. 1 is the radius of the 
 #neighbourhood
R.SetInterpolator(sitk.sitkLinear)
R.SetOptimizerAsGradientDescent(
learningRate=1, numberOfIterations=50, estimateLearningRate=R.EachIteration
) #changed learning rate to 1
R.SetOptimizerScalesFromPhysicalShift()
T = sitk.CenteredTransformInitializer(
new_im_fix,
im_mov,
sitk.AffineTransform(new_im_fix.GetDimension()),
sitk.CenteredTransformInitializerFilter.GEOMETRY,
) #we used 3D. This is where the main problem is perhaps
R.SetInitialTransform(T)


R.SetOptimizerScalesFromPhysicalShift()
#R.SetShrinkFactorsPerLevel(shrinkFactors=[2, 1, 1]) #Changed smooting and shrink factor
#R.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
R.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()


Tx_1pass = R.Execute(new_im_fix, im_mov)

Tx_2pass = R.Execute(new_im_fix, im_mov)
rval_2pass = R.GetMetricValue()

print("  " + R.GetOptimizerStopConditionDescription())
print("  Metric value: {:.3f}".format(rval_2pass))

#composite_transform = R.Execute(sitk.Cast(im_fix, sitk.sitkFloat32), 
                                           #sitk.Cast(im_mov, sitk.sitkFloat32))

registered_image = sitk.Resample(img, img_fix, Tx_2pass) 
reg_array = sitk.GetArrayFromImage(registered_image)