rewrote but test failing...

.vscode/settings.json

@@ -0,0 +1,5 @@
+{
+  "python.testing.pytestArgs": ["tests"],
+  "python.testing.unittestEnabled": false,
+  "python.testing.pytestEnabled": true
+}

examples/stitching_example.py

@@ -12,40 +12,40 @@
 props_file_path = path.join(script_path, "../tests/data/testimages_props.csv")
 images = np.load(image_file_path)
 props = pd.read_csv(props_file_path, index_col=0)
-rows = props["col"].to_list()
-cols = props["row"].to_list()
+rows = props["row"].to_list()
+cols = props["col"].to_list()
 
 print(images.shape)
 # must be 3-dim, with each dimension meaning (tile_index,x,y)
 print(rows)
-# the row indices for each tile index. for example, [1,1,2,2,2,...]
+# the row indices (direction in the last dimension) for each tile index. for example, [1,1,2,2,2,...]
 print(cols)
-# the column indices for each tile index. for example, [2,3,1,2,3,...]
+# the column indices (direction in the second last dinemsion) for each tile index. for example, [2,3,1,2,3,...]
 
-result_df, _ = m2stitch.stitch_images(images, rows, cols)
+result_df, _ = m2stitch.stitch_images(images, rows, cols, row_col_transpose=False)
+# note: the previous default row_col_transpose=True is deprecated and will be removed
 
-print(result_df["x_pos"])
-# the absolute x positions of the tiles
 print(result_df["y_pos"])
-# the absolute y positions of the tiles
-
+# the absolute y positions (second last dimension) of the tiles
+print(result_df["x_pos"])
+# the absolute x positions (last dimension) of the tiles
 
 # stitching example
-result_df["x_pos2"] = result_df["x_pos"] - result_df["x_pos"].min()
 result_df["y_pos2"] = result_df["y_pos"] - result_df["y_pos"].min()
+result_df["x_pos2"] = result_df["x_pos"] - result_df["x_pos"].min()
 
-size_x = images.shape[1]
-size_y = images.shape[2]
+size_y = images.shape[-2]
+size_x = images.shape[-1]
 
 stitched_image_size = (
-    result_df["x_pos2"].max() + size_x,
     result_df["y_pos2"].max() + size_y,
+    result_df["x_pos2"].max() + size_x,
 )
 stitched_image = np.zeros_like(images, shape=stitched_image_size)
 for i, row in result_df.iterrows():
     stitched_image[
-        row["x_pos2"] : row["x_pos2"] + size_x,
         row["y_pos2"] : row["y_pos2"] + size_y,
+        row["x_pos2"] : row["x_pos2"] + size_x,
     ] = images[i]
 
 result_image_file_path = path.join(script_path, "stitched_image.npy")

src/m2stitch/_translation_computation.py

@@ -37,17 +37,14 @@ def pcm(image1: NumArray, image2: NumArray) -> FloatArray:
 
 
 def multi_peak_max(
-    PCM: FloatArray, n: int = 2
+    PCM: FloatArray,
 ) -> Tuple[IntArray, IntArray, FloatArray]:
     """Find the first to n th largest peaks in PCM.
 
     Parameters
     ---------
     PCM : np.ndarray
         the peak correlation matrix
-    n : Int
-        the number of the peaks
-
 
     Returns
     -------
@@ -58,9 +55,9 @@ def multi_peak_max(
     vals : np.ndarray
         the values of the peaks
     """
-    row, col = np.unravel_index(np.argsort(PCM.ravel()), PCM.shape)
-    vals: FloatArray = PCM[row[-n:][::-1], col[-n:][::-1]]
-    return row[-n:][::-1], col[-n:][::-1], vals
+    ys, xs = np.unravel_index(np.argsort(PCM.ravel()), PCM.shape)
+    vals: FloatArray = PCM[ys[::-1], xs[::-1]]
+    return ys[::-1], xs[::-1], vals
 
 
 def ncc(image1: NumArray, image2: NumArray) -> Float:
@@ -89,7 +86,7 @@ def ncc(image1: NumArray, image2: NumArray) -> Float:
     return n / d
 
 
-def extract_overlap_subregion(image: NumArray, x: Int, y: Int) -> NumArray:
+def extract_overlap_subregion(image: NumArray, _y: Int, _x: Int) -> NumArray:
     """Extract the overlapping subregion of the image.
 
     Parameters
@@ -100,37 +97,57 @@ def extract_overlap_subregion(image: NumArray, x: Int, y: Int) -> NumArray:
         the x position
     y : Int
         the y position
-
     Returns
     -------
     subimage : np.ndarray
         the extracted subimage
     """
     sizeY = image.shape[0]
     sizeX = image.shape[1]
-    assert (np.abs(x) < sizeY) and (np.abs(y) < sizeX)
-    xstart = int(max(0, min(x, sizeY, key=int), key=int))
-    xend = int(max(0, min(x + sizeY, sizeY, key=int), key=int))
-    ystart = int(max(0, min(y, sizeX, key=int), key=int))
-    yend = int(max(0, min(y + sizeX, sizeX, key=int), key=int))
+    assert (np.abs(_y) < sizeY) and (np.abs(_x) < sizeX)
+    # clip x to (0, size_Y)
+    xstart = int(max(0, min(_y, sizeY, key=int), key=int))
+    # clip x+sizeY to (0, size_Y)
+    xend = int(max(0, min(_y + sizeY, sizeY, key=int), key=int))
+    ystart = int(max(0, min(_x, sizeX, key=int), key=int))
+    yend = int(max(0, min(_x + sizeX, sizeX, key=int), key=int))
     return image[xstart:xend, ystart:yend]
 
 
 def interpret_translation(
-    image1: NumArray, image2: npt.NDArray, xin: Int, yin: Int
+    image1: NumArray,
+    image2: npt.NDArray,
+    yins: IntArray,
+    xins: IntArray,
+    y_min: Int,
+    y_max: Int,
+    x_min: Int,
+    x_max: Int,
+    n: int = 2,
 ) -> Tuple[float, int, int]:
     """Interpret the translation to find the translation with heighest ncc.
+    The candidates are ... (xin, sizeX-xin) * (+1,-1)
 
     Parameters
     ---------
     image1 : np.ndarray
         the first image (the dimension must be 2)
     image2 : np.ndarray
         the second image (the dimension must be 2)
-    xin : Int
-        the x position estimated by PCM
-    yin : Int
-        the y position estimated by PCM
+    yins : IntArray
+        the y positions estimated by PCM
+    xins : IntArray
+        the x positions estimated by PCM
+    y_min : Int
+        the minimum arrowed y (second last dim.) position of the peak
+    y_max : Int
+        the maximum arrowed y (second last dim.) position of the peak
+    x_min : Int
+        the minimum arrowed x (last dim.) position of the peak
+    x_max : Int
+        the maximum arrowed x (last dim.) position of the peak
+    n : Int
+        the number of the valid peaks to test
 
     Returns
     -------
@@ -145,20 +162,38 @@ def interpret_translation(
     assert image2.ndim == 2
     assert np.array_equal(image1.shape, image2.shape)
     _ncc = -np.infty
-    x = 0
-    y = 0
+    _y = 0
+    _x = 0
     sizeY = image1.shape[0]
     sizeX = image1.shape[1]
-    assert 0 <= xin and xin < sizeY
-    assert 0 <= yin and yin < sizeX
-    xmags = [xin, sizeY - xin] if xin > 0 else [xin]
-    ymags = [yin, sizeX - yin] if yin > 0 else [yin]
-    for xmag, ymag, xsign, ysign in itertools.product(xmags, ymags, [-1, +1], [-1, +1]):
-        subI1 = extract_overlap_subregion(image1, (xmag * xsign), (ymag * ysign))
-        subI2 = extract_overlap_subregion(image2, -(xmag * xsign), -(ymag * ysign))
-        ncc_val = ncc(subI1, subI2)
-        if ncc_val > _ncc:
-            _ncc = float(ncc_val)
-            x = int(xmag * xsign)
-            y = int(ymag * ysign)
-    return _ncc, x, y
+    peak_counts = 0
+
+    for yin, xin in zip(yins, xins):
+        peak_counted = False
+        assert 0 <= yin and yin < sizeY
+        assert 0 <= xin and xin < sizeX
+        ymags = [yin, sizeY - yin] if yin > 0 else [yin]
+        xmags = [xin, sizeX - xin] if xin > 0 else [xin]
+        for ymag, xmag, ysign, xsign in itertools.product(
+            ymags, xmags, [-1, +1], [-1, +1]
+        ):
+            yval = ymag * ysign
+            xval = xmag * xsign
+            if (y_min <= yval) & (yval <= y_max) & (x_min <= xval) & (xval <= x_max):
+                peak_counted = True
+                subI1 = extract_overlap_subregion(
+                    image1, (ymag * ysign), (xmag * xsign)
+                )
+                subI2 = extract_overlap_subregion(
+                    image2, -(ymag * ysign), -(xmag * xsign)
+                )
+                ncc_val = ncc(subI1, subI2)
+                if ncc_val > _ncc:
+                    _ncc = float(ncc_val)
+                    _y = int(ymag * ysign)
+                    _x = int(xmag * ysign)
+        if peak_counted:
+            peak_counts += 1
+            if peak_counts >= n:
+                break
+    return _ncc, _y, _x

src/m2stitch/stitching.py

@@ -138,12 +138,13 @@ def get_index(col, row):
             grid[f"{dimension}_pos_init_guess"] = position_initial_guess[:, j]
         for direction, dimension in itertools.product(["top", "left"], ["y", "x"]):
             for ind, g in grid.iterrows():
-                if g[direction] is not None:
-                    g2 = grid.loc[g[direction]]
-                    grid.loc[ind, f"{direction}_{dimension}_init_guess"] = (
-                        g[f"{dimension}_pos_init_guess"]
-                        - g2[f"{dimension}_pos_init_guess"]
-                    )
+                i1 = g[direction]
+                if pd.isna(i1):
+                    continue
+                g2 = grid.loc[i1]
+                grid.loc[ind, f"{direction}_{dimension}_init_guess"] = (
+                    g2[f"{dimension}_pos_init_guess"] - g[f"{dimension}_pos_init_guess"]
+                )
 
     ###### translationComputation ######
     for direction in ["top", "left"]:
@@ -155,13 +156,31 @@ def get_index(col, row):
             image2 = images[i2]
 
             PCM = pcm(image1, image2).real
+
+            if position_initial_guess is not None:
+                g2 = grid.loc[i1]
+
+                def get_lims(dimension, size):
+                    val = g[f"{direction}_{dimension}_init_guess"]
+                    r = size * overlap_diff_threshold / 100.0
+                    return (val - r, val + r)
+
+                lims = np.array(
+                    [
+                        get_lims(dimension, size)
+                        for dimension, size in zip("yx", [sizeY, sizeX])
+                    ]
+                )
+            else:
+                lims = np.array([[-sizeY, sizeY], [-sizeX, sizeX]])
             found_peaks = list(zip(*multi_peak_max(PCM)))
 
-            interpreted_peaks = []
-            for r, c, _ in found_peaks:
-                interpreted_peaks.append(interpret_translation(image1, image2, r, c))
-            max_peak = interpreted_peaks[np.argmax(np.array(interpreted_peaks)[:, 0])]
-            for j, key in enumerate(["ncc", "x", "y"]):
+            yins, xins, _ = zip(*found_peaks)
+            max_peak = interpret_translation(
+                image1, image2, yins, xins, *lims[0], *lims[1]
+            )
+            #            max_peak = interpreted_peaks[np.argmax(np.array(interpreted_peaks)[:, 0])]
+            for j, key in enumerate(["ncc", "y", "x"]):
                 grid.loc[i2, f"{direction}_{key}_first"] = max_peak[j]
 
     prob_uniform_n, mu_n, sigma_n = compute_image_overlap(

tests/test_stitching.py

@@ -27,3 +27,23 @@ def test_stitching(test_image_path: Tuple[npt.NDArray, pd.DataFrame]) -> None:
     assert np.array_equal(result_df.index, props.index)
     assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 2
     assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 2
+
+
+def test_stitching_initial_guess(
+    test_image_path: Tuple[npt.NDArray, pd.DataFrame]
+) -> None:
+    testimages, props = test_image_path
+    """It exits with a status code of zero."""
+    rows = props["row"].to_list()
+    cols = props["col"].to_list()
+    initial_guess = np.array(props[["y_pos", "x_pos"]].values)
+    result_df, _ = stitch_images(
+        testimages,
+        rows,
+        cols,
+        position_initial_guess=initial_guess,
+        row_col_transpose=False,
+    )
+    assert np.array_equal(result_df.index, props.index)
+    assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 2
+    assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 2