Alignment of the image stack acquired during FIB tomography is essential for accurate representations of the dataset. Small variations in the instrumentation over time, along with the continuous movement of the sample relative to the electron beam with repeated milling need to be corrected. This article will detail various methods to detect the shift in position between images, called image registration, along with methods to move the images between frames, called image warping.
This article will show methods using the Python language to process an example image dataset. To begin coding, one needs to import the packages used into the Python script. In these examples, OpenCV (cv2) and Scikit-Image (si) will be used for image manipulation, Matplotlib (plt) used for visualization, Numpy (np) for matrix manipulation, and the time package will also be used to show the speed of the various registration methods.
import cv2
import numpy as np
import matplotlib.pyplot as plt
import skimage as si
import scipy as sp
import timeTo show the image registration techniques in code, we’ll use data obtained from the Electron Microscopy Data Bank. The images in question are from Cryo-FIB-SEM data on Chlamydomonas reinhardtii cells obtained by Sven Klumpe. To begin, we’ll load in slices 5 and 6 from the dataset and visualize them.
target = si.io.imread('005.tif')
template = si.io.imread('006.tif')
fig,ax = plt.subplots(1,2)
ax[0].imshow(target, cmap = 'gray')
ax[1].imshow(template, cmap = 'gray')
_ = [x.set_axis_off() for x in ax]
ax[0].set_title('Slice 5')
ax[1].set_title('Slice 6')
Before the different alignment methods are detailed, a quick detail on the various motion types in images is in order. Most of the time in FIB tomography, images are only displaced in the x and y directions due to stage drift. In computer vision, movement in the x and y directions only is known as a translation. If the image rotates in the x,y plane along with a translation, the motion is Euclidean. If the image if sheared, elongating a direction, the motion is considered Affine. Finally, the most complicated image motion type is a perspective shift, also known as a homography. Many of the alignment methods below work only for transformations or Euclidean motion. For more details on motion types, see the tutorials from Scikit-Image and Towards Data Science.
Normalized Cross-Correlation (NCC) #
The most often used alignment method is image cross-correlation for images that only contain a translation. Two adjacent images in a z-stack are translated in the x and y direction relative to each other and the matching pixels in the z direction used to calculate the correlation between the images by various objective functions.. The cross-correlation function’s maximum is the translation where the two images are most similar. Most often the images are normalized to reduce the effect of contrast and brightness differences between the two images. While this algorithm gives accurate results, it also can be slow to run if the images are large. Many of the various implementations of cross-correlation can be found in OpenCV’s template matching function.
To perform a normalized cross-correlation (NCC) calculation, OpenCV’s matchTemplate function is used. The function call takes in three arguments, the image that will be aligned (target), the image to align to (template), and a method for calculating the difference between the images (method). For normalized cross-correlation, the TM_CCORR_NORMED method is used. In these examples, we’ll use align Slice 5 to Slice 6, and refer to the stationary image as the template image and the image to be shifted the target image.
For the full NCC calculation, the template image needs to be padded with the size of the target image using the copyMakeBorder function. The padding is necessary so that the target image can be translated to the full width of the template image size. The function call requires the image to be padded (template), and the size of the padding on the top, bottom, left, and right sides of the image, and finally the border method. The padding method uses the cv2.BORDER_CONSTANT method to add values (zero in this example) as the last function call variable to the edge of the images.
ncc_start = time.time() #Timing the opreation
#Pad template to target can be slide across entire image
template_pad = cv2.copyMakeBorder(template, target.shape[0] , target.shape[0],
target.shape[1], target.shape[1],
cv2.BORDER_CONSTANT, None, 0)
corr = cv2.matchTemplate(template_pad,target,cv2.TM_CCOEFF_NORMED)
min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(corr)
#Calculate the translation
target_translate = np.array(max_loc)[::-1] - np.array(target.shape)
ncc_end = time.time()Next, the padded template image and the target image are passed to the cv2.matchTemplate function, along with our objective function, TM_CCORR_NORMED. The corr variable above contains the correlation coefficients between the two images as a function of the translation of the target image. As shown in Figure 2 below, the correlation coefficient has a large peak where the match between the two images is the greatest. One can find the maximum values in the correlation matrix with the minMaxLoc function, which is visualized as the red dot in Figure 2.
fig2,ax2 = plt.subplots(1,1)
ax2.imshow(corr)
#Annotate the max overlap
ax2.scatter(max_loc[0],max_loc[1], marker = '*', color = 'r', s = 30)
Finally, to calculate the translation, one subtracts out the size of the target image to calculate the relative shift in the images:
#Calculate the translation
target_translate = np.array(max_loc)[::-1] - np.array(target.shape)
#Print the results
print(f'Found Shift with NCC: {target_translate}')
ncc_end = time.time()
print(f'Time for OpenCV NCC is {ncc_end-ncc_start} s')In the example above, the translation was a shift to the right of 49 pixels and a shift down of 52 pixels. The time required to perform the calculation was 11.8 s.
Phase Cross-Correlation (PCC) #
Normalized cross-correlation above look for differences in real space intensity cross-correlation functions. Phase cross-correlation looks for maximum overlap of the images in the frequency domain. The two images are converted to a frequency representation using a Fourier transform and the two frequency domain images are then compared by calculating the cross-power spectrum which is at a maximum when the frequencies match between images. Similar to normalized cross-correlation, PCC works only on image translation motion.
The shift to the frequency domain results in calculations that are faster than normalized cross-correlation with large images. The trade-off to shifting to the frequency domain is that the spatial direction is lost, so normal correlations in all 4 translations (+- in x and y) directions need to be calculated to find the correct shift direction. Phase cross-correlation is implemented in OpenCV’s phaseCorrelate function or with Scikit-Image’s phase_cross_correlation function. The Scikit-Image function now also includes an argument to “disabiguate” the shift direction, providing easier registration.
Here we’ll examine the method manually to show the complete registration process. To begin, both the template image and the the target image are converted to the frequency domain with the Fast Fourier Transform (FFT) with Numpy.
target_fft = np.fft.fft2(target)
template_fft = np.fft.fft2(template)
#Plot the images and the ffts
fig4,ax4 = plt.subplots(2,2)
ax4[0,0].imshow(target, cmap = 'gray')
ax4[0,1].imshow(template, cmap = 'gray')
#FFT Mangitude Spectra
ax4[1,0].imshow(np.log(np.abs(np.fft.fftshift(target_fft))), cmap = 'gray')
ax4[1,1].imshow(np.log(np.abs(np.fft.fftshift(template_fft))), cmap = 'gray')
titles = ['Slice 5', 'Slice 6',
'Slice 5 FFT','Slice 6 FFT']
for n,ax in enumerate(ax4.ravel()):
ax.set_axis_off()
ax.set_title(titles[n], size =14)
along with the corresponding FFT power spectra (bottom row).
To find the area with the best overlap in phase space, the cross-power spectrum is calculated by taking element-wise product of the template image with the complex conjugate of the target image, and normalized to the absolute value of the element-wise products of the two FFT images. Normally, a very small number (epsilon) is added to the results to prevent zero division during the normalization. The maximum of the cross-power spectrum is the magnitude shift of the target image versus the template image. In this example the found shift was [41 53].
#Calculate the cross-power spectrum of the ffts
epsilon = np.finfo(target_fft.real.dtype).eps
cps = abs(np.fft.ifft2(
(template_fft * target_fft.conj()) /
(abs(template_fft) * abs(target_fft) + epsilon)
))
#Find the maximum in the cps, the image shift
cps_max = np.array(np.where(cps == cps.max())).T[0]
#Plot the cross-power spectrum
fig5,ax5 = plt.subplots(1,2, figsize = (8,3))
ax5[0].imshow(cps, cmap = 'viridis', origin='lower')
ax5[1].imshow(cps, cmap = 'viridis', origin='lower')
#Zoom in on the maximum in the right figure
ax5[1].set_ylim(0,300)
ax5[1].set_xlim(0,500)
#Annotate the max
ax5[0].scatter(cps_max[1],cps_max[0], marker = '*', color = 'r', s = 30)
ax5[1].scatter(cps_max[1],cps_max[0], marker = '*', color = 'r', s = 30)
fig5.suptitle('Cross Power Spectra')
fig5.tight_layout()
fig5.suptitle('Cross Power Spectra')
fig5.tight_layout()
Red star indicates the maximum of the cross-power spectrum.
The shifts found with phase cross-correlation only determine the magnitude of the shift, not the direction. Thus, it is necessary to calculate the direction of the shift by calculating the correlation between the images by shifting the target image in the ± x and ± y directions. The correlations can be calculated using Numpy’s corrcoef function, along with Scipy’s ndimage.shift function to move the target image.
#measure the absolute correlation for each of the four shifts
directions = [(-1,-1),(-1,1),(1,-1),(1,1)]
real_shift_coeff = []
for shift in directions*cps_max:
#Shift the target image
targ_shift = sp.ndimage.shift(target,shift, mode = 'grid-wrap')
#Find the correlation coeff
coeff = np.corrcoef(targ_shift.ravel(),template.ravel())
real_shift_coeff.append((shift,coeff[0,1]))
real_shift_coeff = sorted(real_shift_coeff, key = lambda x: x[1],
reverse = True)The translation found from phase cross-correlation was a shift to the right of 41 pixels and a shift down of 53 pixels. The time required to perform the calculation was 20.5 s. Of course, this same set of code can be automatically performed using Scikit-Image’s phase_cross_correlation function in a single call:
shift = si.registration.phase_cross_correlation(template, target,
disambiguate = True)With Scikit-Image, the PCC only took 5.4 s, and found the same shift as the manual PCC method.
Enhanced Cross-correlation #
To remove some of the issues with brightness and contrast variation between images, the ‘enhanced cross-correlation‘ (ECC) function was proposed. The ECC function can be used in OpenCV’s findTransformECC and works best when the x-y offset between images is small. Another feature of ECC is that it is able to account for changes in rotation and tilt between images, thus ECC can be calculated using a variety of image motion models: translation, Euclidean, affine, and homography. Of course, looking at motion methods with more degrees of freedom will require more time for optimization. The code below is adapted from LearnOpenCV’s ECC tutorial.
To begin, we need to setup the parameter needed for the ECC calculation. The warp mode is the OpenCV motion model, which can be selected from the four types mentioned above. Next, a warp matrix is initialized as a starting guess for the translation. The np.eye method returns values of 1 along the diagonal of the matrix, which defaults to no movement between images. The size of the warp matrix depends on the motion mode, and should be either a [2×3] for translations, Euclidean, or affine transformations, and [3×3] for homographies. Finally, termination criteria are set for the end of the minimization with a maximum number of iterations, and a early termination if the change in ECC value between two iterations of the calculation is smaller than the termination_eps.
warp_mode = cv2.MOTION_TRANSLATION # Define the motion model
warp_matrix = np.eye(2, 3, dtype=np.float32) #Define Start warp matrix
number_of_iterations = 500 # Specify the number of iterations.
# Specify the threshold of the increment
# in the correlation coefficient between two iterations
termination_eps = 1e-10;
# Define termination criteria
criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT,
number_of_iterations, termination_eps)Running the ECC algorithm uses a call to the cv2.findTransformECC method, and passing in the template image, target image, warp matrix, warp mode, and termination criteria. The function then returns the found correlation coefficient and the optimized warp matrix.
# Run the ECC algorithm. The results are stored in warp_matrix.
(cc, warp_matrix) = cv2.findTransformECC(template,target,warp_matrix,
warp_mode, criteria)
To apply the warp transformation the cv2.warpAffine is used, unless the motion model is homography which uses cv2.warpPerspective. The call the the warp function takes in the image to be warped, the warp matrix that is applied, the desired shape of the the output image. Additionally, flags are passed that describe how on how to interpolate pixels, and the WARP_INVERSE_MAP which tells the function that the transformation matrix will be applied to shift the image onto the source image.
# Use warpAffine for Translation, Euclidean and Affine
sz = im1.shape # Find size of image1
im2_aligned = cv2.warpAffine(im2, warp_matrix, (sz[1],sz[0]),
flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)#Pic of the transformation
#Visualize the results
fig2,ax2 = plt.subplots(2,2)
ax2[0,0].imshow(im1, cmap = 'gray')
ax2[0,1].imshow(im2_aligned, cmap = 'gray')
diff_unalign = cv2.subtract(im1,im2)
diff_align = cv2.subtract(im1, im2_aligned)
ax2[1,0].imshow(diff_unalign, vmin = 0, vmax = 255, cmap = 'gray')
ax2[1,1].imshow(diff_align, vmin = 0, vmax = 255, cmap = 'gray')
titles = ['Slice 1', 'Slice 2 Aligned to Slice 1',
'Difference between Slice 1 and Slice 2',
'Difference between Slice 1 and Aligned Slice 2']
for n,ax in enumerate(ax2.ravel()):
ax.set_axis_off()
ax.set_title(titles[n])As stated above, ECC works best when the translation between images is small and fails completely if the motion is large. One method to overcome this deficiency is to downscale the images to lower resolution and find the ECC warp matrix on the lower resolution image. Using the lower resolution warp matrix as the input for the findTransformECC method in the high resolution image allows the algorithm to start nearest the maximum of the correlation field and can result in finding the best warp matrix.
Feature Matching #
While not used frequently in FIB tomography due to the typical high overlap of images between stack slices, feature matching and alignment can be useful under some situations. If there are significant distortions between images due to lens changes or sample changes or other imaging artifacts, normal cross-correlation methods may fail to provide good alignment. In that case, feature matching may be able to align the images, but takes significantly more time to compute. To learn more about feature matching, see the tutorial at OpenCV.
