Use difference of gaussians for alignment

lens-biophotonics · May 2, 2022 · f4f919f · f4f919f
1 parent aad74ee
commit f4f919f
Show file tree

Hide file tree

Showing 4 changed files with 53 additions and 197 deletions.
diff --git a/scripts/align_debug.py b/scripts/align_debug.py
@@ -6,10 +6,10 @@
 import numpy as np
 import tifffile as tiff
 
-import cv2 as cv
-
 from zetastitcher.io.inputfile import InputFile
 
+from zetastitcher.align.dog import align_dog
+
 
 def to_dtype(x, dtype):
     x = np.rint(x) if np.issubdtype(dtype, np.integer) else x
@@ -76,19 +76,16 @@ def stitch(aname, bname, z_frame, axis, overlap, max_shift_z=20,
         bframe = np.rot90(bframe, axes=(-1, -2))
     b_roi = bframe[..., :overlap - max_shift_y, :].astype(np.float32)
 
-    padding = [(0, 0), (0, 0), (max_shift_x, max_shift_x)]
-
-    a_roi = np.pad(a_roi, padding, 'constant').astype(np.float32)
-
     tiff.imsave('aslice.tiff', a_roi.astype(np.float32))
     tiff.imsave('bframe.tiff', b_roi.astype(np.float32))
 
-    output_shape = np.array(a_roi.shape) - np.array(b_roi.shape) + 1
+    output_shape = np.array(a_roi.shape) + np.array((0, 0, 2 * max_shift_x)) - np.array(b_roi.shape) + 1
     output_shape[0] = a_roi.shape[0]
     xcorr = np.zeros(output_shape)
 
     for i in range(xcorr.shape[0]):
-        xcorr[i] = cv.matchTemplate(a_roi[i], b_roi[0], cv.TM_CCOEFF_NORMED)
+        cc, max_loc = align_dog(a_roi[i], b_roi[0], 0, max_shift_x)
+        xcorr[i] = cc
     tiff.imsave('xcorr.tiff', xcorr.astype(np.float32))
 
     shift = list(np.unravel_index(np.argmax(xcorr), xcorr.shape))

diff --git a/zetastitcher/align/__main__.py b/zetastitcher/align/__main__.py
@@ -19,9 +19,9 @@
 from zetastitcher.io.inputfile import InputFile
 from zetastitcher.align.filematrix import FileMatrix
 from zetastitcher.align.xcorr_filematrix import XcorrFileMatrix
+from zetastitcher.align.dog import align_dog
 from zetastitcher.fuse import absolute_positions
 from zetastitcher.fuse.__main__ import ABS_MODE_MAXIMUM_SCORE
-from .normxcorr import normxcorr2_cv
 
 from zetastitcher.version import __version__
 
@@ -168,9 +168,6 @@ def worker(item, overlap_dict, channel, max_dz, max_dy, max_dx):
         aslice = np.rot90(aslice, axes=(-1, -2))
     aslice = aslice[..., -(overlap + max_dy):, :]
 
-    padding = [(0, 0), (0, 0), (max_dx, max_dx)]
-    aslice = np.pad(aslice, padding, 'constant')
-
     bframe = b.zslice_idx(z_frame, copy=True)
     if axis == 2:
         bframe = np.rot90(bframe, axes=(-1, -2))
@@ -179,7 +176,14 @@ def worker(item, overlap_dict, channel, max_dz, max_dy, max_dx):
     aslice = aslice.astype(np.float32)
     bframe = bframe.astype(np.float32)
 
-    xcorr = normxcorr2_cv(aslice, bframe)
+    output_shape = np.array(aslice.shape) + np.array((0, 0, 2 * max_dx)) - np.array(bframe.shape) + 1
+    output_shape[0] = aslice.shape[0]
+    xcorr = np.zeros(output_shape)
+
+    for i in range(xcorr.shape[0]):
+        cc, max_loc = align_dog(aslice[i], bframe[0], 0, max_dx)
+        xcorr[i] = cc
+
     shift = list(np.unravel_index(np.argmax(xcorr), xcorr.shape))
     score = xcorr[tuple(shift)]
     if score < 0 or score > 1:

diff --git a/zetastitcher/align/dog.py b/zetastitcher/align/dog.py
@@ -0,0 +1,39 @@
+import cv2 as cv
+
+import numpy as np
+
+
+def twoD_gaussian_kernel(ksize, sigma):
+    filter = cv.getGaussianKernel(ksize, sigma)
+    return cv.mulTransposed(filter, False)
+
+
+def dog(image, ksize=100, sigma1=5, sigma2=10):
+    filter1 = twoD_gaussian_kernel(ksize, sigma1)
+    filter2 = twoD_gaussian_kernel(ksize, sigma2)
+    cv.normalize(filter1, filter1, 1, 0, cv.NORM_L1)
+    cv.normalize(filter2, filter2, 1, 0, cv.NORM_L1)
+    filter = filter1 - filter2
+
+    return cv.filter2D(image, -1, filter)
+
+
+def crossCorr(image1, image2, padding_y, padding_x):
+    temp1 = image1.astype(np.float32)
+    temp2 = image2.astype(np.float32)
+
+    padding = (padding_y, padding_y, padding_x, padding_x)
+
+    padded = cv.copyMakeBorder(temp1, *padding, cv.BORDER_CONSTANT)
+
+    return cv.matchTemplate(padded, temp2, cv.TM_CCORR_NORMED)
+
+
+def align_dog(i1, i2, padding_y, padding_x):
+    dog1 = dog(i1)
+    dog2 = dog(i2)
+
+    cc = crossCorr(dog1, dog2, padding_y, padding_x)
+
+    min_val, max_val, min_loc, max_loc = cv.minMaxLoc(cc)
+    return cc, max_loc
diff --git a/zetastitcher/align/normxcorr.py b/zetastitcher/align/normxcorr.py