Search code examples
pythonimage-processingnumpydicomimage-registration

Easy way to calculate the shift between two arrays? Python


I have two dicom images and can compute the mean squared error of the images using the code below. However, there can be an inherent shift in one image compared to another (if my imager is slightly misaligned). Is there an easy way to compute a shift of two numpy arrays?

I have tried shifting the array by a few pixels each direction and computing the minimum MSQ. However, this is not robust enough. Any help or advice would be very much appreciated!

import numpy as np
import dicom

#first image
ds = dicom.read_file("U:\\temp\\1.dcm")
array1 = ds.pixel_array
#second image
ds1 = dicom.read_file("U:\\temp\\5.dcm")
array2 = ds1.pixel_array

#shifting image by a number of pixels in any direction
arr1 = np.roll(array2, 100,  axis=1)
imshow(arr1)

def mse(imageA, imageB):
    # the 'Mean Squared Error' between the two images is the
    # sum of the squared difference between the two images;
    err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2)
    err /= float(imageA.shape[0] * imageA.shape[1])

    # return the MSE, the lower the error, the more "similar"
    # the two images are
    return err

first_try = mse(array1, array2)
second_try = mse(arr1, array2)

Solution

  • If you are sure that the images will be exactly identical except for this shift, think the best solution is scipy.signal.correlate2d:

    import numpy as np
    from scipy.signal import correlate2d
    
    a = np.random.random((200, 200))
    b = np.roll(np.roll(a, 15, axis=0),-33, axis=1)
    
    corr = correlate2d(a, b)
    
    shift = np.where(corr==corr.max())
    shift = (shift[0]%a.shape[0], shift[1]%a.shape[1])
    

    This gives the correct values:

    (array([184]), array([32]))