From 0696f24f18a9994476e775bd0cbaf21461079e63 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Thu, 10 Feb 2022 22:54:09 +0900 Subject: [PATCH 01/31] added args and refactored row cols --- examples/stitching_example.py | 4 ++-- src/m2stitch/__main__.py | 4 ++-- src/m2stitch/_stage_model.py | 10 +++++----- src/m2stitch/stitching.py | 28 +++++++++++++++++++++------- tests/test_stitching.py | 4 ++-- 5 files changed, 32 insertions(+), 18 deletions(-) diff --git a/examples/stitching_example.py b/examples/stitching_example.py index 953033d..dddf020 100755 --- a/examples/stitching_example.py +++ b/examples/stitching_example.py @@ -12,8 +12,8 @@ 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["row"].to_list() -cols = props["col"].to_list() +rows = props["col"].to_list() +cols = props["row"].to_list() print(images.shape) # must be 3-dim, with each dimension meaning (tile_index,x,y) diff --git a/src/m2stitch/__main__.py b/src/m2stitch/__main__.py index bc9c96f..ac0696f 100644 --- a/src/m2stitch/__main__.py +++ b/src/m2stitch/__main__.py @@ -12,8 +12,8 @@ def main() -> None: # testimages, props = test_image_path # """It exits with a status code of zero.""" -# rows = props["row"].to_list() -# cols = props["col"].to_list() +# rows = props["col"].to_list() +# cols = props["row"].to_list() # result_df, _ = stitch_images(testimages, rows, cols) diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index 238f1a9..f633f43 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -223,7 +223,7 @@ def compute_repeatability( if np.any(grid["top_valid2"]): xs = grid[grid["top_valid2"]]["top_x_first"] rx_top = np.ceil((xs.max() - xs.min()) / 2) - _, yss = zip(*grid[grid["top_valid2"]].groupby("col")["top_y_first"]) + _, yss = zip(*grid[grid["top_valid2"]].groupby("row")["top_y_first"]) ry_top = np.ceil(np.max([np.max(ys) - np.min(ys) for ys in yss]) / 2) r_top = max(rx_top, ry_top) else: @@ -232,7 +232,7 @@ def compute_repeatability( if np.any(grid["left_valid2"]): ys = grid[grid["left_valid2"]]["left_y_first"] ry_left = np.ceil((ys.max() - ys.min()) / 2) - _, xss = zip(*grid[grid["left_valid2"]].groupby("row")["left_x_first"]) + _, xss = zip(*grid[grid["left_valid2"]].groupby("col")["left_x_first"]) rx_left = np.ceil(np.max([np.max(xs) - np.min(xs) for xs in xss]) / 2) r_left = max(ry_left, rx_left) else: @@ -256,7 +256,7 @@ def filter_by_repeatability(grid: pd.DataFrame, r: Float) -> pd.DataFrame: grid : pd.DataFrame the updated dataframe for the grid position """ - for _, grp in grid.groupby("row"): + for _, grp in grid.groupby("col"): isvalid = grp["top_valid2"].astype(bool) if not any(isvalid): grid.loc[grp.index, "top_valid3"] = False @@ -268,7 +268,7 @@ def filter_by_repeatability(grid: pd.DataFrame, r: Float) -> pd.DataFrame: & grp["top_y_first"].between(medy - r, medy + r) & (grp["top_ncc_first"] > 0.5) ) - for _, grp in grid.groupby("col"): + for _, grp in grid.groupby("row"): isvalid = grp["left_valid2"] if not any(isvalid): grid.loc[grp.index, "left_valid3"] = False @@ -303,7 +303,7 @@ def replace_invalid_translations(grid: pd.DataFrame) -> pd.DataFrame: grid.loc[isvalid, f"{direction}_{key}_second"] = grid.loc[ isvalid, f"{direction}_{key}_first" ] - for direction, rowcol in zip(["top", "left"], ["row", "col"]): + for direction, rowcol in zip(["top", "left"], ["col", "row"]): for _, grp in grid.groupby(rowcol): isvalid = grp[f"{direction}_valid3"].astype(bool) if any(isvalid): diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index 0e6fd1d..d67869a 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -29,6 +29,8 @@ def stitch_images( rows: Optional[Sequence[Any]] = None, cols: Optional[Sequence[Any]] = None, position_indices: Optional[NumArray] = None, + position_initial_guess: Optional[NumArray] = None, + overlap_diff_threshold: Float = 10, pou: Float = 3, overlap_prob_uniform_threshold: Float = 80, full_output: bool = False, @@ -52,6 +54,13 @@ def stitch_images( the dimensions corresponds to (image, index) ignored if rows and cols are not None. + position_initial_guess : np.ndarray, optional + the initial guess for the positions of the images, in the unit of pixels. + + overlap_diff_threshold : 10 + the allowed difference from the initial guess, in percentage of the image size. + ignored if position_initial_guess is None + pou : Float, default 3 the "percent overlap uncertainty" parameter @@ -89,21 +98,26 @@ def stitch_images( position_indices = np.array(position_indices) assert images.shape[0] == position_indices.shape[0] assert position_indices.shape[1] == images.ndim - 1 + if position_initial_guess is not None: + position_initial_guess = np.array(position_initial_guess) + assert images.shape[0] == position_indices.shape[0] + assert position_initial_guess.shape[1] == images.ndim - 1 assert 0 <= overlap_prob_uniform_threshold and overlap_prob_uniform_threshold <= 100 - _rows, _cols = position_indices.T + assert 0 <= overlap_diff_threshold and overlap_diff_threshold <= 100 + _cols, _rows = position_indices.T W, H = images.shape[1:] grid = pd.DataFrame( { - "row": _rows, "col": _cols, + "row": _rows, }, - index=np.arange(len(_rows)), + index=np.arange(len(_cols)), ) def get_index(row, col): - df = grid[(grid["row"] == row) & (grid["col"] == col)] + df = grid[(grid["col"] == row) & (grid["row"] == col)] assert len(df) < 2 if len(df) == 1: return df.index[0] @@ -111,10 +125,10 @@ def get_index(row, col): return None grid["top"] = grid.apply( - lambda g: get_index(g["row"] - 1, g["col"]), axis=1 + lambda g: get_index(g["col"] - 1, g["row"]), axis=1 ).astype(pd.Int32Dtype()) grid["left"] = grid.apply( - lambda g: get_index(g["row"], g["col"] - 1), axis=1 + lambda g: get_index(g["col"], g["row"] - 1), axis=1 ).astype(pd.Int32Dtype()) ###### translationComputation ###### @@ -177,4 +191,4 @@ def get_index(row, col): if full_output: return grid, prop_dict else: - return grid[["col", "row", "y_pos", "x_pos"]], prop_dict + return grid[["row", "col", "y_pos", "x_pos"]], prop_dict diff --git a/tests/test_stitching.py b/tests/test_stitching.py index a7d0359..970e5e0 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -21,8 +21,8 @@ def test_image_path(shared_datadir: str) -> Tuple[npt.NDArray, pd.DataFrame]: def test_stitching(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() + rows = props["col"].to_list() + cols = props["row"].to_list() result_df, _ = stitch_images(testimages, rows, cols) assert np.array_equal(result_df.index, props.index) assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 2 From 621d9a2b3c3a48cf3f5bf0e9ae3478eaa2f07fa4 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Thu, 10 Feb 2022 23:08:18 +0900 Subject: [PATCH 02/31] added some initial guess method --- src/m2stitch/stitching.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index d67869a..4a31a47 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -1,4 +1,5 @@ """This module provides microscope image stitching with the algorithm by MIST.""" +import itertools import warnings from typing import Any from typing import Optional @@ -116,8 +117,8 @@ def stitch_images( index=np.arange(len(_cols)), ) - def get_index(row, col): - df = grid[(grid["col"] == row) & (grid["row"] == col)] + def get_index(col, row): + df = grid[(grid["col"] == col) & (grid["row"] == row)] assert len(df) < 2 if len(df) == 1: return df.index[0] @@ -131,6 +132,19 @@ def get_index(row, col): lambda g: get_index(g["col"], g["row"] - 1), axis=1 ).astype(pd.Int32Dtype()) + ### dimension order ... m.y.x + if position_initial_guess is not None: + for j, dimension in enumerate(["y", "x"]): + 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"] + ) + ###### translationComputation ###### for direction in ["top", "left"]: for i2, g in tqdm(grid.iterrows(), total=len(grid)): From 96ed6e9c6dd140dfc996f4311a1c3ba356628143 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Thu, 10 Feb 2022 23:12:41 +0900 Subject: [PATCH 03/31] H and W refactoring --- LICENSE.rst | 2 +- src/m2stitch/_constrained_refinement.py | 12 ++++++++--- src/m2stitch/_stage_model.py | 27 ++++++++++++++---------- src/m2stitch/_translation_computation.py | 26 +++++++++++------------ src/m2stitch/stitching.py | 16 ++++++++------ 5 files changed, 49 insertions(+), 34 deletions(-) diff --git a/LICENSE.rst b/LICENSE.rst index 24abdc4..745fe95 100644 --- a/LICENSE.rst +++ b/LICENSE.rst @@ -29,6 +29,6 @@ The original license for MIST (https://github.com/usnistgov/MIST/blob/27a8787/LI NIST-developed software is provided by NIST as a public service. You may use, copy and distribute copies of the software in any medium, provided that you keep intact this entire notice. You may improve, modify and create derivative works of the software or any portion of the software, and you may copy and distribute such modifications or works. Modified works should carry a notice stating that you changed the software and should note the date and nature of any such change. Please explicitly acknowledge the National Institute of Standards and Technology as the source of the software. -NIST-developed software is expressly provided "AS IS." NIST MAKES NO WARRANTY OF ANY KIND, EXPRESS, IMPLIED, IN FACT OR ARISING BY OPERATION OF LAW, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NON-INFRINGEMENT AND DATA ACCURACY. NIST NEITHER REPRESENTS NOR WARRANTS THAT THE OPERATION OF THE SOFTWARE WILL BE UNINTERRUPTED OR ERROR-FREE, OR THAT ANY DEFECTS WILL BE CORRECTED. NIST DOES NOT WARRANT OR MAKE ANY REPRESENTATIONS REGARDING THE USE OF THE SOFTWARE OR THE RESULTS THEREOF, INCLUDING BUT NOT LIMITED TO THE CORRECTNESS, ACCURACY, RELIABILITY, OR USEFULNESS OF THE SOFTWARE. +NIST-developed software is expressly provided "AS IS." NIST MAKES NO WARRANTY OF ANY KIND, EXPRESS, IMPLIED, IN FACT OR ARISING BY OPERATION OF LAsizeY, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NON-INFRINGEMENT AND DATA ACCURACY. NIST NEITHER REPRESENTS NOR WARRANTS THAT THE OPERATION OF THE SOFTWARE WILL BE UNINTERRUPTED OR ERROR-FREE, OR THAT ANY DEFECTS WILL BE CORRECTED. NIST DOES NOT WARRANT OR MAKE ANY REPRESENTATIONS REGARDING THE USE OF THE SOFTWARE OR THE RESULTS THEREOF, INCLUDING BUT NOT LIMITED TO THE CORRECTNESS, ACCURACY, RELIABILITY, OR USEFULNESS OF THE SOFTWARE. You are solely responsible for determining the appropriateness of using and distributing the software and you assume all risks associated with its use, including but not limited to the risks and costs of program errors, compliance with applicable laws, damage to or loss of data, programs or equipment, and the unavailability or interruption of operation. This software is not intended to be used in any situation where a failure could cause risk of injury or damage to property. The software developed by NIST employees is not subject to copyright protection within the United States. diff --git a/src/m2stitch/_constrained_refinement.py b/src/m2stitch/_constrained_refinement.py index 1f6e72c..591b3de 100644 --- a/src/m2stitch/_constrained_refinement.py +++ b/src/m2stitch/_constrained_refinement.py @@ -90,7 +90,7 @@ def refine_translations(images: NumArray, grid: pd.DataFrame, r: Float) -> pd.Da continue image1 = images[i1] image2 = images[i2] - W, H = image1.shape + sizeY, sizeX = image1.shape def overlap_ncc(params): x, y = params @@ -103,8 +103,14 @@ def overlap_ncc(params): int(g[f"{direction}_y_second"]), ] limits = [ - [max(-W + 1, init_values[0] - r), min(W - 1, init_values[0] + r)], - [max(-H + 1, init_values[1] - r), min(H - 1, init_values[1] + r)], + [ + max(-sizeY + 1, init_values[0] - r), + min(sizeY - 1, init_values[0] + r), + ], + [ + max(-sizeX + 1, init_values[1] - r), + min(sizeX - 1, init_values[1] + r), + ], ] values, ncc_value = find_local_max_integer_constrained( overlap_ncc, np.array(init_values), np.array(limits) diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index f633f43..dbafb43 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -63,8 +63,8 @@ def compute_inv_liklihood( def compute_image_overlap( grid: pd.DataFrame, direction: str, - W: Int, - H: Int, + sizeY: Int, + sizeX: Int, max_stall_count: Int = 100, prob_uniform_threshold: Float = 80, ) -> Tuple[float, ...]: @@ -76,9 +76,9 @@ def compute_image_overlap( the dataframe for the grid position, with columns "{top|left}_{x|y}_first" direction : str the direction of the overlap, either of "top" or "left" - W : Int + sizeY : Int the image width - H : Int + sizeX : Int the image height max_stall_count : Int, optional the maximum count of optimization, by default 100 @@ -96,9 +96,9 @@ def compute_image_overlap( when direction is not in ["top","left"], raises ValueError """ if direction == "top": - T = grid["top_y_first"].values / H * 100 + T = grid["top_y_first"].values / sizeX * 100 elif direction == "left": - T = grid["left_x_first"].values / W * 100 + T = grid["left_x_first"].values / sizeY * 100 else: raise ValueError("direction must be either of top or left") T = T[np.isfinite(T)] @@ -185,7 +185,12 @@ def filter_outliers(T: pd.Series, isvalid: pd.Series, w: Float = 1.5) -> pd.Seri def compute_repeatability( - grid: pd.DataFrame, overlap_n: Float, overlap_w: Float, W: Int, H: Int, pou: Float + grid: pd.DataFrame, + overlap_n: Float, + overlap_w: Float, + sizeY: Int, + sizeX: Int, + pou: Float, ) -> Tuple[pd.DataFrame, Float]: """Compute the repeatability of the stage motion. @@ -197,9 +202,9 @@ def compute_repeatability( the estimated overlap for top direction overlap_w : Float the estimated overlap for left direction - W : Int + sizeY : Int the width of the height images - H : Int + sizeX : Int the width of the tile images pou : Float the percentile margin for error, by default 3 @@ -212,10 +217,10 @@ def compute_repeatability( the repeatability of the stage motion """ grid["top_valid1"] = filter_by_overlap_and_correlation( - grid["top_y_first"], grid["top_ncc_first"], overlap_n, H, pou + grid["top_y_first"], grid["top_ncc_first"], overlap_n, sizeX, pou ) grid["left_valid1"] = filter_by_overlap_and_correlation( - grid["left_x_first"], grid["left_ncc_first"], overlap_w, W, pou + grid["left_x_first"], grid["left_ncc_first"], overlap_w, sizeY, pou ) grid["top_valid2"] = filter_outliers(grid["top_y_first"], grid["top_valid1"]) grid["left_valid2"] = filter_outliers(grid["left_x_first"], grid["left_valid1"]) diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index 1fed241..6271fae 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -106,13 +106,13 @@ def extract_overlap_subregion(image: NumArray, x: Int, y: Int) -> NumArray: subimage : np.ndarray the extracted subimage """ - W = image.shape[0] - H = image.shape[1] - assert (np.abs(x) < W) and (np.abs(y) < H) - xstart = int(max(0, min(x, W, key=int), key=int)) - xend = int(max(0, min(x + W, W, key=int), key=int)) - ystart = int(max(0, min(y, H, key=int), key=int)) - yend = int(max(0, min(y + H, H, key=int), key=int)) + 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)) return image[xstart:xend, ystart:yend] @@ -147,12 +147,12 @@ def interpret_translation( _ncc = -np.infty x = 0 y = 0 - W = image1.shape[0] - H = image1.shape[1] - assert 0 <= xin and xin < W - assert 0 <= yin and yin < H - xmags = [xin, W - xin] if xin > 0 else [xin] - ymags = [yin, H - yin] if yin > 0 else [yin] + 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)) diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index 4a31a47..7a1f3bc 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -107,7 +107,7 @@ def stitch_images( assert 0 <= overlap_diff_threshold and overlap_diff_threshold <= 100 _cols, _rows = position_indices.T - W, H = images.shape[1:] + sizeY, sizeX = images.shape[1:] grid = pd.DataFrame( { @@ -165,18 +165,22 @@ def get_index(col, row): grid.loc[i2, f"{direction}_{key}_first"] = max_peak[j] prob_uniform_n, mu_n, sigma_n = compute_image_overlap( - grid, "top", W, H, prob_uniform_threshold=overlap_prob_uniform_threshold + grid, "top", sizeY, sizeX, prob_uniform_threshold=overlap_prob_uniform_threshold ) overlap_n = 100 - mu_n prob_uniform_w, mu_w, sigma_w = compute_image_overlap( - grid, "left", W, H, prob_uniform_threshold=overlap_prob_uniform_threshold + grid, + "left", + sizeY, + sizeX, + prob_uniform_threshold=overlap_prob_uniform_threshold, ) overlap_w = 100 - mu_w overlap_n = np.clip(overlap_n, pou, 100 - pou) overlap_w = np.clip(overlap_w, pou, 100 - pou) - grid, r = compute_repeatability(grid, overlap_n, overlap_w, W, H, pou) + grid, r = compute_repeatability(grid, overlap_n, overlap_w, sizeY, sizeX, pou) grid = filter_by_repeatability(grid, r) grid = replace_invalid_translations(grid) @@ -186,8 +190,8 @@ def get_index(col, row): grid = compute_final_position(grid, tree) prop_dict = { - "W": W, - "H": H, + "W": sizeY, + "H": sizeX, "overlap_top": overlap_n, "overlap_left": overlap_w, "overlap_top_results": { From 5e235ce8afe608b9d5565340ee6650e2f40867a1 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Fri, 11 Feb 2022 07:39:13 +0900 Subject: [PATCH 04/31] test passing --- tests/data/testimages_props.csv | 2 +- tests/test_stitching.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/data/testimages_props.csv b/tests/data/testimages_props.csv index e09ef1b..220fbc0 100644 --- a/tests/data/testimages_props.csv +++ b/tests/data/testimages_props.csv @@ -1,4 +1,4 @@ -,row,col,x_pos,y_pos +,col,row,x_pos,y_pos 0,2,0,1,3330 1,3,0,0,5000 2,3,1,1349,5001 diff --git a/tests/test_stitching.py b/tests/test_stitching.py index 970e5e0..bfa4840 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -21,9 +21,9 @@ def test_image_path(shared_datadir: str) -> Tuple[npt.NDArray, pd.DataFrame]: def test_stitching(test_image_path: Tuple[npt.NDArray, pd.DataFrame]) -> None: testimages, props = test_image_path """It exits with a status code of zero.""" - rows = props["col"].to_list() - cols = props["row"].to_list() - result_df, _ = stitch_images(testimages, rows, cols) + rows = props["row"].to_list() + cols = props["col"].to_list() + result_df, _ = stitch_images(testimages, rows, cols, 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 From 946488c13e7f3a79f8132d7b5b8da541edfe72a8 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sat, 12 Feb 2022 15:13:19 +0900 Subject: [PATCH 05/31] rewrote but test failing... --- .vscode/settings.json | 5 ++ examples/stitching_example.py | 28 +++--- src/m2stitch/_translation_computation.py | 103 +++++++++++++++-------- src/m2stitch/stitching.py | 41 ++++++--- tests/test_stitching.py | 20 +++++ 5 files changed, 138 insertions(+), 59 deletions(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..d969f96 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "python.testing.pytestArgs": ["tests"], + "python.testing.unittestEnabled": false, + "python.testing.pytestEnabled": true +} diff --git a/examples/stitching_example.py b/examples/stitching_example.py index dddf020..b1645b8 100755 --- a/examples/stitching_example.py +++ b/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") diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index 6271fae..39ada5d 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -37,7 +37,7 @@ 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. @@ -45,9 +45,6 @@ def multi_peak_max( --------- 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,7 +97,6 @@ def extract_overlap_subregion(image: NumArray, x: Int, y: Int) -> NumArray: the x position y : Int the y position - Returns ------- subimage : np.ndarray @@ -108,18 +104,29 @@ def extract_overlap_subregion(image: NumArray, x: Int, y: Int) -> NumArray: """ 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 --------- @@ -127,10 +134,20 @@ def interpret_translation( 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 diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index 7a1f3bc..3d411ca 100644 --- a/src/m2stitch/stitching.py +++ b/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( diff --git a/tests/test_stitching.py b/tests/test_stitching.py index bfa4840..79bf0dd 100644 --- a/tests/test_stitching.py +++ b/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 From fdd2f156f1a258e68213aab3ad172d10f150e914 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sat, 12 Feb 2022 19:55:50 +0900 Subject: [PATCH 06/31] Revert "rewrote but test failing..." This reverts commit 946488c13e7f3a79f8132d7b5b8da541edfe72a8. --- .vscode/settings.json | 5 -- examples/stitching_example.py | 28 +++--- src/m2stitch/_translation_computation.py | 103 ++++++++--------------- src/m2stitch/stitching.py | 41 +++------ tests/test_stitching.py | 20 ----- 5 files changed, 59 insertions(+), 138 deletions(-) delete mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index d969f96..0000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,5 +0,0 @@ -{ - "python.testing.pytestArgs": ["tests"], - "python.testing.unittestEnabled": false, - "python.testing.pytestEnabled": true -} diff --git a/examples/stitching_example.py b/examples/stitching_example.py index b1645b8..dddf020 100755 --- a/examples/stitching_example.py +++ b/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["row"].to_list() -cols = props["col"].to_list() +rows = props["col"].to_list() +cols = props["row"].to_list() print(images.shape) # must be 3-dim, with each dimension meaning (tile_index,x,y) print(rows) -# the row indices (direction in the last dimension) for each tile index. for example, [1,1,2,2,2,...] +# the row indices for each tile index. for example, [1,1,2,2,2,...] print(cols) -# the column indices (direction in the second last dinemsion) for each tile index. for example, [2,3,1,2,3,...] +# the column indices for each tile index. for example, [2,3,1,2,3,...] -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 +result_df, _ = m2stitch.stitch_images(images, rows, cols) -print(result_df["y_pos"]) -# the absolute y positions (second last dimension) of the tiles print(result_df["x_pos"]) -# the absolute x positions (last dimension) of the tiles +# the absolute x positions of the tiles +print(result_df["y_pos"]) +# the absolute y positions of the tiles + # stitching example -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() +result_df["y_pos2"] = result_df["y_pos"] - result_df["y_pos"].min() -size_y = images.shape[-2] -size_x = images.shape[-1] +size_x = images.shape[1] +size_y = images.shape[2] stitched_image_size = ( - result_df["y_pos2"].max() + size_y, result_df["x_pos2"].max() + size_x, + result_df["y_pos2"].max() + size_y, ) stitched_image = np.zeros_like(images, shape=stitched_image_size) for i, row in result_df.iterrows(): stitched_image[ - row["y_pos2"] : row["y_pos2"] + size_y, row["x_pos2"] : row["x_pos2"] + size_x, + row["y_pos2"] : row["y_pos2"] + size_y, ] = images[i] result_image_file_path = path.join(script_path, "stitched_image.npy") diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index 39ada5d..6271fae 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -37,7 +37,7 @@ def pcm(image1: NumArray, image2: NumArray) -> FloatArray: def multi_peak_max( - PCM: FloatArray, + PCM: FloatArray, n: int = 2 ) -> Tuple[IntArray, IntArray, FloatArray]: """Find the first to n th largest peaks in PCM. @@ -45,6 +45,9 @@ def multi_peak_max( --------- PCM : np.ndarray the peak correlation matrix + n : Int + the number of the peaks + Returns ------- @@ -55,9 +58,9 @@ def multi_peak_max( vals : np.ndarray the values of the peaks """ - ys, xs = np.unravel_index(np.argsort(PCM.ravel()), PCM.shape) - vals: FloatArray = PCM[ys[::-1], xs[::-1]] - return ys[::-1], xs[::-1], vals + 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 def ncc(image1: NumArray, image2: NumArray) -> Float: @@ -86,7 +89,7 @@ def ncc(image1: NumArray, image2: NumArray) -> Float: return n / d -def extract_overlap_subregion(image: NumArray, _y: Int, _x: Int) -> NumArray: +def extract_overlap_subregion(image: NumArray, x: Int, y: Int) -> NumArray: """Extract the overlapping subregion of the image. Parameters @@ -97,6 +100,7 @@ def extract_overlap_subregion(image: NumArray, _y: Int, _x: Int) -> NumArray: the x position y : Int the y position + Returns ------- subimage : np.ndarray @@ -104,29 +108,18 @@ def extract_overlap_subregion(image: NumArray, _y: Int, _x: Int) -> NumArray: """ sizeY = image.shape[0] sizeX = image.shape[1] - 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)) + 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)) return image[xstart:xend, ystart:yend] def interpret_translation( - image1: NumArray, - image2: npt.NDArray, - yins: IntArray, - xins: IntArray, - y_min: Int, - y_max: Int, - x_min: Int, - x_max: Int, - n: int = 2, + image1: NumArray, image2: npt.NDArray, xin: Int, yin: Int ) -> Tuple[float, int, int]: """Interpret the translation to find the translation with heighest ncc. - The candidates are ... (xin, sizeX-xin) * (+1,-1) Parameters --------- @@ -134,20 +127,10 @@ def interpret_translation( the first image (the dimension must be 2) image2 : np.ndarray the second image (the dimension must be 2) - 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 + xin : Int + the x position estimated by PCM + yin : Int + the y position estimated by PCM Returns ------- @@ -162,38 +145,20 @@ def interpret_translation( assert image2.ndim == 2 assert np.array_equal(image1.shape, image2.shape) _ncc = -np.infty - _y = 0 - _x = 0 + x = 0 + y = 0 sizeY = image1.shape[0] sizeX = image1.shape[1] - 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 + 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 diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index 3d411ca..7a1f3bc 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -138,13 +138,12 @@ 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(): - 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"] - ) + 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"] + ) ###### translationComputation ###### for direction in ["top", "left"]: @@ -156,31 +155,13 @@ 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))) - 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"]): + 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"]): grid.loc[i2, f"{direction}_{key}_first"] = max_peak[j] prob_uniform_n, mu_n, sigma_n = compute_image_overlap( diff --git a/tests/test_stitching.py b/tests/test_stitching.py index 79bf0dd..bfa4840 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -27,23 +27,3 @@ 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 From 27e1a1da60d07c23b9cc512762e78de902f4bf76 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sat, 12 Feb 2022 19:56:31 +0900 Subject: [PATCH 07/31] added pytest --- .vscode/settings.json | 1 + 1 file changed, 1 insertion(+) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..0967ef4 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1 @@ +{} From 21097572697c7a2f5ae095b1c4229fc7283e07e7 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sat, 12 Feb 2022 19:58:31 +0900 Subject: [PATCH 08/31] settings updated --- .vscode/settings.json | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.vscode/settings.json b/.vscode/settings.json index 0967ef4..d969f96 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1 +1,5 @@ -{} +{ + "python.testing.pytestArgs": ["tests"], + "python.testing.unittestEnabled": false, + "python.testing.pytestEnabled": true +} From 74b772df9f4abf958c7335cd062b11970351c994 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sat, 12 Feb 2022 20:15:25 +0900 Subject: [PATCH 09/31] refacting translation computation --- src/m2stitch/_translation_computation.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index 6271fae..a8ffa0e 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -89,18 +89,17 @@ 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 --------- image : np.ndarray the image (the dimension must be 2) - x : Int - the x position y : Int - the y position - + the y (second last dim.) position + x : Int + the x (last dim.) position Returns ------- subimage : np.ndarray @@ -108,11 +107,13 @@ def extract_overlap_subregion(image: NumArray, x: Int, y: Int) -> NumArray: """ 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] From 424575540acac9d5f10fddf608486454c82d954d Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sat, 12 Feb 2022 22:41:37 +0900 Subject: [PATCH 10/31] test working --- src/m2stitch/_translation_computation.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index a8ffa0e..b381c8d 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -118,7 +118,7 @@ def extract_overlap_subregion(image: NumArray, y: Int, x: Int) -> NumArray: def interpret_translation( - image1: NumArray, image2: npt.NDArray, xin: Int, yin: Int + image1: NumArray, image2: npt.NDArray, yin: Int, xin: Int ) -> Tuple[float, int, int]: """Interpret the translation to find the translation with heighest ncc. @@ -150,16 +150,18 @@ def interpret_translation( y = 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)) + 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 + subI1 = extract_overlap_subregion(image1, yval, xval) + subI2 = extract_overlap_subregion(image2, -yval, -xval) 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 + y = int(yval) + x = int(xval) + return _ncc, y, x From 5a500fbecd17750805b96ac609e72d73e253fd2e Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sat, 12 Feb 2022 22:47:59 +0900 Subject: [PATCH 11/31] interpret_translation refactor --- src/m2stitch/_translation_computation.py | 49 ++++++++++++++---------- src/m2stitch/stitching.py | 8 +--- 2 files changed, 31 insertions(+), 26 deletions(-) diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index b381c8d..e746ae2 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -118,7 +118,7 @@ def extract_overlap_subregion(image: NumArray, y: Int, x: Int) -> NumArray: def interpret_translation( - image1: NumArray, image2: npt.NDArray, yin: Int, xin: Int + image1: NumArray, image2: npt.NDArray, yins: IntArray, xins: IntArray, n: Int = 2 ) -> Tuple[float, int, int]: """Interpret the translation to find the translation with heighest ncc. @@ -128,10 +128,12 @@ def interpret_translation( 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 + xins : IntArray + the x positions estimated by PCM + yins : IntArray + the y positions estimated by PCM + n : Int + the number of peaks to check, default is 2. Returns ------- @@ -150,18 +152,25 @@ def interpret_translation( y = 0 sizeY = image1.shape[0] sizeX = image1.shape[1] - 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 - subI1 = extract_overlap_subregion(image1, yval, xval) - subI2 = extract_overlap_subregion(image2, -yval, -xval) - ncc_val = ncc(subI1, subI2) - if ncc_val > _ncc: - _ncc = float(ncc_val) - y = int(yval) - x = int(xval) - return _ncc, y, x + checked_num = 0 + for yin, xin in zip(yins, xins): + 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 + subI1 = extract_overlap_subregion(image1, yval, xval) + subI2 = extract_overlap_subregion(image2, -yval, -xval) + ncc_val = ncc(subI1, subI2) + if ncc_val > _ncc: + _ncc = float(ncc_val) + y = int(yval) + x = int(xval) + checked_num += 1 + if checked_num >= n: + break + return _ncc, y, x diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index 7a1f3bc..a4d976c 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -155,12 +155,8 @@ def get_index(col, row): image2 = images[i2] PCM = pcm(image1, image2).real - 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])] + yins, xins, _ = multi_peak_max(PCM) + max_peak = interpret_translation(image1, image2, yins, xins) for j, key in enumerate(["ncc", "x", "y"]): grid.loc[i2, f"{direction}_{key}_first"] = max_peak[j] From fd8afa1e1429eb944c57d54e3955404c858ea51a Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sat, 12 Feb 2022 22:51:57 +0900 Subject: [PATCH 12/31] test working... --- src/m2stitch/_translation_computation.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index e746ae2..01221e8 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -36,18 +36,13 @@ def pcm(image1: NumArray, image2: NumArray) -> FloatArray: return np.fft.ifft2(FC / np.abs(FC)).real.astype(np.float32) -def multi_peak_max( - PCM: FloatArray, n: int = 2 -) -> Tuple[IntArray, IntArray, FloatArray]: +def multi_peak_max(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 ------- @@ -59,8 +54,8 @@ def multi_peak_max( 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 + vals: FloatArray = PCM[row[::-1], col[::-1]] + return row[::-1], col[::-1], vals def ncc(image1: NumArray, image2: NumArray) -> Float: @@ -173,4 +168,4 @@ def interpret_translation( checked_num += 1 if checked_num >= n: break - return _ncc, y, x + return _ncc, y, x From e1ca8d1d82c0d3b47aed721e048664e36f990d42 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 08:57:29 +0900 Subject: [PATCH 13/31] test working... --- src/m2stitch/_translation_computation.py | 42 +++++++++++++++++------- src/m2stitch/stitching.py | 32 ++++++++++++++---- 2 files changed, 56 insertions(+), 18 deletions(-) diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index 01221e8..efb2516 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -113,7 +113,15 @@ def extract_overlap_subregion(image: NumArray, y: Int, x: Int) -> NumArray: def interpret_translation( - image1: NumArray, image2: npt.NDArray, yins: IntArray, xins: IntArray, n: Int = 2 + image1: NumArray, + image2: npt.NDArray, + yins: IntArray, + xins: IntArray, + ymin: Int, + ymax: Int, + xmin: Int, + xmax: Int, + n: Int = 2, ) -> Tuple[float, int, int]: """Interpret the translation to find the translation with heighest ncc. @@ -123,10 +131,18 @@ def interpret_translation( the first image (the dimension must be 2) image2 : np.ndarray the second image (the dimension must be 2) - xins : IntArray - the x positions estimated by PCM yins : IntArray the y positions estimated by PCM + xins : IntArray + the x positions estimated by PCM + ymin : Int + the minimum value of y (second last dim.) + ymax : Int + the maximum value of y (second last dim.) + xmin : Int + the minimum value of x (last dim.) + xmax : Int + the maximum value of x (last dim.) n : Int the number of peaks to check, default is 2. @@ -149,6 +165,7 @@ def interpret_translation( sizeX = image1.shape[1] checked_num = 0 for yin, xin in zip(yins, xins): + ischecked = False assert 0 <= yin and yin < sizeY assert 0 <= xin and xin < sizeX ymags = [yin, sizeY - yin] if yin > 0 else [yin] @@ -158,14 +175,17 @@ def interpret_translation( ): yval = ymag * ysign xval = xmag * xsign - subI1 = extract_overlap_subregion(image1, yval, xval) - subI2 = extract_overlap_subregion(image2, -yval, -xval) - ncc_val = ncc(subI1, subI2) - if ncc_val > _ncc: - _ncc = float(ncc_val) - y = int(yval) - x = int(xval) - checked_num += 1 + if (ymin <= yval) and (yval <= ymax) and (xmin <= xval) and (xval <= xmax): + ischecked = True + subI1 = extract_overlap_subregion(image1, yval, xval) + subI2 = extract_overlap_subregion(image2, -yval, -xval) + ncc_val = ncc(subI1, subI2) + if ncc_val > _ncc: + _ncc = float(ncc_val) + y = int(yval) + x = int(xval) + if ischecked: + checked_num += 1 if checked_num >= n: break return _ncc, y, x diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index a4d976c..6c86348 100644 --- a/src/m2stitch/stitching.py +++ b/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"] = ( + g[f"{dimension}_pos_init_guess"] - g2[f"{dimension}_pos_init_guess"] + ) ###### translationComputation ###### for direction in ["top", "left"]: @@ -155,8 +156,25 @@ def get_index(col, row): image2 = images[i2] PCM = pcm(image1, image2).real + if position_initial_guess is not None: + + 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]]) yins, xins, _ = multi_peak_max(PCM) - max_peak = interpret_translation(image1, image2, yins, xins) + max_peak = interpret_translation( + image1, image2, yins, xins, *lims[0], *lims[1] + ) for j, key in enumerate(["ncc", "x", "y"]): grid.loc[i2, f"{direction}_{key}_first"] = max_peak[j] From 4dfc80b42e71d84341bcb918c90c985f1de1fb9e Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 11:58:24 +0900 Subject: [PATCH 14/31] test working --- src/m2stitch/_translation_computation.py | 49 +++++++++++++++--------- tests/test_stitching.py | 25 +++++++++++- 2 files changed, 53 insertions(+), 21 deletions(-) diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index efb2516..162aaa8 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -158,25 +158,40 @@ def interpret_translation( assert image1.ndim == 2 assert image2.ndim == 2 assert np.array_equal(image1.shape, image2.shape) + sizeY = image1.shape[0] + sizeX = image1.shape[1] + assert np.all(0 <= yins) and np.all(yins < sizeY) + assert np.all(0 <= xins) and np.all(xins < sizeX) + _ncc = -np.infty x = 0 y = 0 - sizeY = image1.shape[0] - sizeX = image1.shape[1] - checked_num = 0 - for yin, xin in zip(yins, xins): - ischecked = 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 + + ymagss = [yins, sizeY - yins] + ymagss[1][ymagss[0] == 0] = 0 + xmagss = [xins, sizeX - xins] + xmagss[1][xmagss[0] == 0] = 0 + + # concatenate all the candidates + poss = [] + for ymags, xmags, ysign, xsign in itertools.product( + ymagss, xmagss, [-1, +1], [-1, +1] + ): + yvals = ymags * ysign + xvals = xmags * xsign + poss.append([yvals, xvals]) + poss = np.array(poss) + valid_ind = ( + (ymin <= poss[:, 0, :]) + & (poss[:, 0, :] <= ymax) + & (xmin <= poss[:, 1, :]) + & (poss[:, 1, :] <= xmax) + ) + assert np.any(valid_ind) + valid_ind = np.any(valid_ind, axis=0) + for pos in np.moveaxis(poss[:, :, valid_ind], -1, 0)[:n]: + for yval, xval in pos: if (ymin <= yval) and (yval <= ymax) and (xmin <= xval) and (xval <= xmax): - ischecked = True subI1 = extract_overlap_subregion(image1, yval, xval) subI2 = extract_overlap_subregion(image2, -yval, -xval) ncc_val = ncc(subI1, subI2) @@ -184,8 +199,4 @@ def interpret_translation( _ncc = float(ncc_val) y = int(yval) x = int(xval) - if ischecked: - checked_num += 1 - if checked_num >= n: - break return _ncc, y, x diff --git a/tests/test_stitching.py b/tests/test_stitching.py index bfa4840..cf13d58 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -21,9 +21,30 @@ def test_image_path(shared_datadir: str) -> Tuple[npt.NDArray, pd.DataFrame]: def test_stitching(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() + rows = props["row"].to_list() result_df, _ = stitch_images(testimages, rows, cols, 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 + assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 2 + + +def test_stitching_init_guess( + test_image_path: Tuple[npt.NDArray, pd.DataFrame] +) -> None: + np.random.seed(0) + testimages, props = test_image_path + """It exits with a status code of zero.""" + cols = props["col"].to_list() + rows = props["row"].to_list() + pos_guess = props[["x_pos", "y_pos"]].values + result_df, _ = stitch_images( + testimages, + rows, + cols, + row_col_transpose=False, + position_initial_guess=pos_guess, + ) + assert np.array_equal(result_df.index, props.index) + assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) <= 3 + assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) <= 3 From 3d2450e17c11e66c13edf000edfe7329990d423e Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 12:08:10 +0900 Subject: [PATCH 15/31] changed _first variables --- src/m2stitch/_stage_model.py | 40 ++++++++++++++++++------------------ src/m2stitch/stitching.py | 2 +- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index dbafb43..3167470 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -96,9 +96,9 @@ def compute_image_overlap( when direction is not in ["top","left"], raises ValueError """ if direction == "top": - T = grid["top_y_first"].values / sizeX * 100 + T = grid["top_x_first"].values / sizeX * 100 elif direction == "left": - T = grid["left_x_first"].values / sizeY * 100 + T = grid["left_y_first"].values / sizeY * 100 else: raise ValueError("direction must be either of top or left") T = T[np.isfinite(T)] @@ -217,27 +217,27 @@ def compute_repeatability( the repeatability of the stage motion """ grid["top_valid1"] = filter_by_overlap_and_correlation( - grid["top_y_first"], grid["top_ncc_first"], overlap_n, sizeX, pou + grid["top_x_first"], grid["top_ncc_first"], overlap_n, sizeX, pou ) grid["left_valid1"] = filter_by_overlap_and_correlation( - grid["left_x_first"], grid["left_ncc_first"], overlap_w, sizeY, pou + grid["left_y_first"], grid["left_ncc_first"], overlap_w, sizeY, pou ) - grid["top_valid2"] = filter_outliers(grid["top_y_first"], grid["top_valid1"]) - grid["left_valid2"] = filter_outliers(grid["left_x_first"], grid["left_valid1"]) + grid["top_valid2"] = filter_outliers(grid["top_x_first"], grid["top_valid1"]) + grid["left_valid2"] = filter_outliers(grid["left_y_first"], grid["left_valid1"]) if np.any(grid["top_valid2"]): - xs = grid[grid["top_valid2"]]["top_x_first"] + xs = grid[grid["top_valid2"]]["top_y_first"] rx_top = np.ceil((xs.max() - xs.min()) / 2) - _, yss = zip(*grid[grid["top_valid2"]].groupby("row")["top_y_first"]) + _, yss = zip(*grid[grid["top_valid2"]].groupby("row")["top_x_first"]) ry_top = np.ceil(np.max([np.max(ys) - np.min(ys) for ys in yss]) / 2) r_top = max(rx_top, ry_top) else: r_top = 0 # better than failing if np.any(grid["left_valid2"]): - ys = grid[grid["left_valid2"]]["left_y_first"] + ys = grid[grid["left_valid2"]]["left_x_first"] ry_left = np.ceil((ys.max() - ys.min()) / 2) - _, xss = zip(*grid[grid["left_valid2"]].groupby("col")["left_x_first"]) + _, xss = zip(*grid[grid["left_valid2"]].groupby("col")["left_y_first"]) rx_left = np.ceil(np.max([np.max(xs) - np.min(xs) for xs in xss]) / 2) r_left = max(ry_left, rx_left) else: @@ -266,11 +266,11 @@ def filter_by_repeatability(grid: pd.DataFrame, r: Float) -> pd.DataFrame: if not any(isvalid): grid.loc[grp.index, "top_valid3"] = False else: - medx = grp[isvalid]["top_x_first"].median() - medy = grp[isvalid]["top_y_first"].median() + medx = grp[isvalid]["top_y_first"].median() + medy = grp[isvalid]["top_x_first"].median() grid.loc[grp.index, "top_valid3"] = ( - grp["top_x_first"].between(medx - r, medx + r) - & grp["top_y_first"].between(medy - r, medy + r) + grp["top_y_first"].between(medx - r, medx + r) + & grp["top_x_first"].between(medy - r, medy + r) & (grp["top_ncc_first"] > 0.5) ) for _, grp in grid.groupby("row"): @@ -278,11 +278,11 @@ def filter_by_repeatability(grid: pd.DataFrame, r: Float) -> pd.DataFrame: if not any(isvalid): grid.loc[grp.index, "left_valid3"] = False else: - medx = grp[isvalid]["left_x_first"].median() - medy = grp[isvalid]["left_y_first"].median() + medx = grp[isvalid]["left_y_first"].median() + medy = grp[isvalid]["left_x_first"].median() grid.loc[grp.index, "left_valid3"] = ( - grp["left_x_first"].between(medx - r, medx + r) - & grp["left_y_first"].between(medy - r, medy + r) + grp["left_y_first"].between(medx - r, medx + r) + & grp["left_x_first"].between(medy - r, medy + r) & (grp["left_ncc_first"] > 0.5) ) return grid @@ -319,10 +319,10 @@ def replace_invalid_translations(grid: pd.DataFrame) -> pd.DataFrame: pd.isna(grid.loc[grp.index[~isvalid], f"{direction}_y_second"]) ) grid.loc[grp.index[~isvalid], f"{direction}_x_second"] = grp[isvalid][ - f"{direction}_x_first" + f"{direction}_y_first" ].median() grid.loc[grp.index[~isvalid], f"{direction}_y_second"] = grp[isvalid][ - f"{direction}_y_first" + f"{direction}_x_first" ].median() grid.loc[grp.index[~isvalid], f"{direction}_ncc_second"] = -1 for direction, xy in itertools.product(["top", "left"], ["x", "y"]): diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index 6c86348..f64f50a 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -175,7 +175,7 @@ def get_lims(dimension, size): max_peak = interpret_translation( image1, image2, yins, xins, *lims[0], *lims[1] ) - for j, key in enumerate(["ncc", "x", "y"]): + 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( From fb13c6866be751cafb2f2e5789f54f4a3c7898a1 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 12:11:45 +0900 Subject: [PATCH 16/31] further refacting --- src/m2stitch/_constrained_refinement.py | 2 +- src/m2stitch/_stage_model.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/m2stitch/_constrained_refinement.py b/src/m2stitch/_constrained_refinement.py index 591b3de..e053db6 100644 --- a/src/m2stitch/_constrained_refinement.py +++ b/src/m2stitch/_constrained_refinement.py @@ -99,8 +99,8 @@ def overlap_ncc(params): return ncc(subI1, subI2) init_values = [ - int(g[f"{direction}_x_second"]), int(g[f"{direction}_y_second"]), + int(g[f"{direction}_x_second"]), ] limits = [ [ diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index 3167470..46cc5db 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -313,15 +313,15 @@ def replace_invalid_translations(grid: pd.DataFrame) -> pd.DataFrame: isvalid = grp[f"{direction}_valid3"].astype(bool) if any(isvalid): assert all( - pd.isna(grid.loc[grp.index[~isvalid], f"{direction}_x_second"]) + pd.isna(grid.loc[grp.index[~isvalid], f"{direction}_y_second"]) ) assert all( - pd.isna(grid.loc[grp.index[~isvalid], f"{direction}_y_second"]) + pd.isna(grid.loc[grp.index[~isvalid], f"{direction}_x_second"]) ) - grid.loc[grp.index[~isvalid], f"{direction}_x_second"] = grp[isvalid][ + grid.loc[grp.index[~isvalid], f"{direction}_y_second"] = grp[isvalid][ f"{direction}_y_first" ].median() - grid.loc[grp.index[~isvalid], f"{direction}_y_second"] = grp[isvalid][ + grid.loc[grp.index[~isvalid], f"{direction}_x_second"] = grp[isvalid][ f"{direction}_x_first" ].median() grid.loc[grp.index[~isvalid], f"{direction}_ncc_second"] = -1 From a96e9e05d03f6c996467e996beb332f903fe5274 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 12:17:44 +0900 Subject: [PATCH 17/31] test working (refactoring done) --- src/m2stitch/_constrained_refinement.py | 8 ++++---- src/m2stitch/_global_optimization.py | 21 ++++++++++----------- tests/test_stitching.py | 4 ++-- 3 files changed, 16 insertions(+), 17 deletions(-) diff --git a/src/m2stitch/_constrained_refinement.py b/src/m2stitch/_constrained_refinement.py index e053db6..7c0f6b0 100644 --- a/src/m2stitch/_constrained_refinement.py +++ b/src/m2stitch/_constrained_refinement.py @@ -115,11 +115,11 @@ def overlap_ncc(params): values, ncc_value = find_local_max_integer_constrained( overlap_ncc, np.array(init_values), np.array(limits) ) - grid.loc[i2, f"{direction}_x"] = values[0] - grid.loc[i2, f"{direction}_y"] = values[1] + grid.loc[i2, f"{direction}_y"] = values[0] + grid.loc[i2, f"{direction}_x"] = values[1] grid.loc[i2, f"{direction}_ncc"] = ncc_value for direction in ["top", "left"]: - for xy in ["x", "y"]: - key = f"{direction}_{xy}" + for dim in "yx": + key = f"{direction}_{dim}" grid[key] = grid[key].astype(pd.Int32Dtype()) return grid diff --git a/src/m2stitch/_global_optimization.py b/src/m2stitch/_global_optimization.py index 0b38ae4..5e482cc 100644 --- a/src/m2stitch/_global_optimization.py +++ b/src/m2stitch/_global_optimization.py @@ -33,8 +33,8 @@ def compute_maximum_spanning_tree(grid: pd.DataFrame) -> nx.Graph: direction=direction, f=i, t=g[direction], - x=g[f"{direction}_x"], y=g[f"{direction}_y"], + x=g[f"{direction}_x"], ) return nx.maximum_spanning_tree(connection_graph) @@ -58,8 +58,8 @@ def compute_final_position( grid : pd.DataFrame the result dataframe for the grid position, with columns "{x|y}_pos" """ - grid.loc[source_index, "x_pos"] = 0 grid.loc[source_index, "y_pos"] = 0 + grid.loc[source_index, "x_pos"] = 0 nodes = [source_index] walked_nodes = [] @@ -72,20 +72,19 @@ def compute_final_position( props["t"] == node ) & (props["f"] == adj) nodes.append(adj) - x_pos = grid.loc[node, "x_pos"] y_pos = grid.loc[node, "y_pos"] + x_pos = grid.loc[node, "x_pos"] if node == props["t"]: - grid.loc[adj, "x_pos"] = x_pos + props["x"] grid.loc[adj, "y_pos"] = y_pos + props["y"] + grid.loc[adj, "x_pos"] = x_pos + props["x"] else: - grid.loc[adj, "x_pos"] = x_pos - props["x"] grid.loc[adj, "y_pos"] = y_pos - props["y"] - assert not any(pd.isna(grid["x_pos"])) - assert not any(pd.isna(grid["y_pos"])) - grid["x_pos"] = grid["x_pos"] - grid["x_pos"].min() - grid["y_pos"] = grid["y_pos"] - grid["y_pos"].min() - grid["x_pos"] = grid["x_pos"].astype(np.int32) - grid["y_pos"] = grid["y_pos"].astype(np.int32) + grid.loc[adj, "x_pos"] = x_pos - props["x"] + for dim in "yx": + k = f"{dim}_pos" + assert not any(pd.isna(grid[k])) + grid[k] = grid[k] - grid[k].min() + grid[k] = grid[k].astype(np.int32) return grid diff --git a/tests/test_stitching.py b/tests/test_stitching.py index cf13d58..181c60b 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -25,8 +25,8 @@ def test_stitching(test_image_path: Tuple[npt.NDArray, pd.DataFrame]) -> None: rows = props["row"].to_list() result_df, _ = stitch_images(testimages, rows, cols, row_col_transpose=False) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 2 - assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 2 + assert np.max(np.abs(result_df["y_pos"] - props["x_pos"])) < 2 + assert np.max(np.abs(result_df["x_pos"] - props["y_pos"])) < 2 def test_stitching_init_guess( From a37375a4d48b4e04f31151e8f52dc435405239b0 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 12:21:11 +0900 Subject: [PATCH 18/31] further refactoring position --- src/m2stitch/stitching.py | 2 ++ tests/data/testimages_props.csv | 2 +- tests/test_stitching.py | 6 +++--- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index f64f50a..d490a97 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -220,6 +220,8 @@ def get_lims(dimension, size): }, "repeatability": r, } + if row_col_transpose: + grid = grid.rename(columns={"x_pos": "y_pos", "y_pos": "x_pos"}) if full_output: return grid, prop_dict else: diff --git a/tests/data/testimages_props.csv b/tests/data/testimages_props.csv index 220fbc0..fb30292 100644 --- a/tests/data/testimages_props.csv +++ b/tests/data/testimages_props.csv @@ -1,4 +1,4 @@ -,col,row,x_pos,y_pos +,col,row,y_pos,x_pos 0,2,0,1,3330 1,3,0,0,5000 2,3,1,1349,5001 diff --git a/tests/test_stitching.py b/tests/test_stitching.py index 181c60b..5b1f6a4 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -25,8 +25,8 @@ def test_stitching(test_image_path: Tuple[npt.NDArray, pd.DataFrame]) -> None: rows = props["row"].to_list() result_df, _ = stitch_images(testimages, rows, cols, row_col_transpose=False) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["x_pos"])) < 2 - assert np.max(np.abs(result_df["x_pos"] - props["y_pos"])) < 2 + assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 2 + assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 2 def test_stitching_init_guess( @@ -37,7 +37,7 @@ def test_stitching_init_guess( """It exits with a status code of zero.""" cols = props["col"].to_list() rows = props["row"].to_list() - pos_guess = props[["x_pos", "y_pos"]].values + pos_guess = props[["y_pos", "x_pos"]].values result_df, _ = stitch_images( testimages, rows, From 4763273161454763fa66e99680753bc7c7681411 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 12:35:31 +0900 Subject: [PATCH 19/31] solving nox problems in safety and mypy --- noxfile.py | 5 ++++- src/m2stitch/_translation_computation.py | 6 +++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/noxfile.py b/noxfile.py index ac4f132..c9432a8 100644 --- a/noxfile.py +++ b/noxfile.py @@ -101,7 +101,10 @@ def safety(session: Session) -> None: """Scan dependencies for insecure packages.""" requirements = session.poetry.export_requirements() session.install("safety") - session.run("safety", "check", "--full-report", f"--file={requirements}") + session.run( + "safety", "check", "--full-report", f"--file={requirements}", "--ignore=44715" + ) + # ignore numpy update for a while as resolved in 1.22 (2022.2.13) @session(python=python_versions) diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index 162aaa8..64222f7 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -173,14 +173,14 @@ def interpret_translation( xmagss[1][xmagss[0] == 0] = 0 # concatenate all the candidates - poss = [] + _poss = [] for ymags, xmags, ysign, xsign in itertools.product( ymagss, xmagss, [-1, +1], [-1, +1] ): yvals = ymags * ysign xvals = xmags * xsign - poss.append([yvals, xvals]) - poss = np.array(poss) + _poss.append([yvals, xvals]) + poss = np.array(_poss) valid_ind = ( (ymin <= poss[:, 0, :]) & (poss[:, 0, :] <= ymax) From 8d4ba346009efe167c8b378bdf89bd42eeda5a59 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 12:37:45 +0900 Subject: [PATCH 20/31] further --- src/m2stitch/_translation_computation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index 64222f7..5572314 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -189,7 +189,7 @@ def interpret_translation( ) assert np.any(valid_ind) valid_ind = np.any(valid_ind, axis=0) - for pos in np.moveaxis(poss[:, :, valid_ind], -1, 0)[:n]: + for pos in np.moveaxis(poss[:, :, valid_ind], -1, 0)[: int(n)]: for yval, xval in pos: if (ymin <= yval) and (yval <= ymax) and (xmin <= xval) and (xval <= xmax): subI1 = extract_overlap_subregion(image1, yval, xval) From dbc678d70e4617fe1e121092db33d5552fc93554 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 13:35:54 +0900 Subject: [PATCH 21/31] solved typeguard problem --- src/m2stitch/stitching.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index d490a97..fa8d38a 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -161,7 +161,7 @@ def get_index(col, row): def get_lims(dimension, size): val = g[f"{direction}_{dimension}_init_guess"] r = size * overlap_diff_threshold / 100.0 - return (val - r, val + r) + return np.round([val - r, val + r]).astype(np.int64) lims = np.array( [ @@ -178,6 +178,8 @@ def get_lims(dimension, size): for j, key in enumerate(["ncc", "y", "x"]): grid.loc[i2, f"{direction}_{key}_first"] = max_peak[j] + # top_displacement= + # left_displacement= prob_uniform_n, mu_n, sigma_n = compute_image_overlap( grid, "top", sizeY, sizeX, prob_uniform_threshold=overlap_prob_uniform_threshold ) From 683cf40897891a2227939a9a1f56324d433a93d3 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:11:31 +0900 Subject: [PATCH 22/31] refactor filtering logic to isolationforest --- src/m2stitch/_constrained_refinement.py | 8 +- src/m2stitch/_global_optimization.py | 4 +- src/m2stitch/_stage_model.py | 218 ++++------------------- src/m2stitch/_translation_computation.py | 2 +- src/m2stitch/stitching.py | 70 ++++---- tests/test_stitching.py | 8 +- 6 files changed, 85 insertions(+), 225 deletions(-) diff --git a/src/m2stitch/_constrained_refinement.py b/src/m2stitch/_constrained_refinement.py index 7c0f6b0..7f2cd03 100644 --- a/src/m2stitch/_constrained_refinement.py +++ b/src/m2stitch/_constrained_refinement.py @@ -74,16 +74,16 @@ def refine_translations(images: NumArray, grid: pd.DataFrame, r: Float) -> pd.Da images : NumArray the tile images grid : pd.DataFrame - the dataframe for the grid position, with columns "{top|left}_{x|y}_second" + the dataframe for the grid position, with columns "{left|top}_{x|y}_second" r : Float the repeatability Returns ------- grid : pd.DataFrame - the refined grid position, with columns "{top|left}_{x|y|ncc}" + the refined grid position, with columns "{left|top}_{x|y|ncc}" """ - for direction in ["top", "left"]: + for direction in ["left", "top"]: for i2, g in tqdm(grid.iterrows(), total=len(grid)): i1 = g[direction] if pd.isna(i1): @@ -118,7 +118,7 @@ def overlap_ncc(params): grid.loc[i2, f"{direction}_y"] = values[0] grid.loc[i2, f"{direction}_x"] = values[1] grid.loc[i2, f"{direction}_ncc"] = ncc_value - for direction in ["top", "left"]: + for direction in ["left", "top"]: for dim in "yx": key = f"{direction}_{dim}" grid[key] = grid[key].astype(pd.Int32Dtype()) diff --git a/src/m2stitch/_global_optimization.py b/src/m2stitch/_global_optimization.py index 5e482cc..5f560b9 100644 --- a/src/m2stitch/_global_optimization.py +++ b/src/m2stitch/_global_optimization.py @@ -12,7 +12,7 @@ def compute_maximum_spanning_tree(grid: pd.DataFrame) -> nx.Graph: ---------- grid : pd.DataFrame the dataframe for the grid position, - with columns "{top|left}_{x|y|ncc|valid3}" + with columns "{left|top}_{x|y|ncc|valid3}" Returns ------- @@ -21,7 +21,7 @@ def compute_maximum_spanning_tree(grid: pd.DataFrame) -> nx.Graph: """ connection_graph = nx.Graph() for i, g in grid.iterrows(): - for direction in ["top", "left"]: + for direction in ["left", "top"]: if not pd.isna(g[direction]): weight = g[f"{direction}_ncc"] if g[f"{direction}_valid3"]: diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index 46cc5db..a48784c 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -1,135 +1,53 @@ import itertools from typing import Tuple -from typing import Union import numpy as np import pandas as pd -from scipy.optimize import minimize +from sklearn.ensemble import IsolationForest from ._typing_utils import Float -from ._typing_utils import FloatArray from ._typing_utils import Int -def calc_liklihood(prob_uniform: Float, mu: Float, sigma: Float, t: Float) -> Float: - """Calculate the liklihood of the translation. - - Parameters - ---------- - prob_uniform : Float - the probability of the uniform error - mu : Float - the mean of the Gaussian distribution - sigma : Float - the stdev of the Gaussian distribution - t : Float - the translation amplitude - - Returns - ------- - liklihood : Float - the calculated liklihood - """ - t2 = -((t - mu) ** 2) / (2 * sigma ** 2) - norm_liklihood = 1.0 / (np.sqrt(2 * np.pi) * sigma) * np.exp(t2) - uniform_liklihood = 1 / 100.0 - p = prob_uniform / 100.0 - return p * uniform_liklihood + (1 - p) * norm_liklihood - - -def compute_inv_liklihood( - params: Union[Tuple[Float, Float, Float], FloatArray], T: FloatArray -) -> Float: - """Compute the negated sum of log liklihood. - - Parameters - ---------- - params : Union[Tuple[Float, Float, Float], FloatArray] - the parameters : prob_uniform, mu, sigma - T : FloatArray - the value of translations - - Returns - ------- - inv_log_likelihood : Float - the negated sum of log liklihood - """ - prob_uniform, mu, sigma = params - return -np.sum( - [np.log(np.abs(calc_liklihood(prob_uniform, mu, sigma, t))) for t in T] - ) - - -def compute_image_overlap( +def compute_image_overlap2( grid: pd.DataFrame, direction: str, sizeY: Int, sizeX: Int, - max_stall_count: Int = 100, - prob_uniform_threshold: Float = 80, -) -> Tuple[float, ...]: +) -> Tuple[Float, Float]: """Compute the value of the image overlap. Parameters ---------- grid : pd.DataFrame - the dataframe for the grid position, with columns "{top|left}_{x|y}_first" + the dataframe for the grid position, with columns "{left|top}_{x|y}_first" direction : str - the direction of the overlap, either of "top" or "left" + the direction of the overlap, either of "left" or "top" sizeY : Int the image width sizeX : Int the image height - max_stall_count : Int, optional - the maximum count of optimization, by default 100 - prob_uniform_threshold : Float, optional - the threshold for the probabillity for uniform error, by default 80 Returns ------- - x : Tuple[float, float, float] - the computed parameters : prob_uniform, mu, sigma + x : Tuple[float, float] + the computed y and x displacement Raises ------ ValueError - when direction is not in ["top","left"], raises ValueError + when direction is not in ["left","top"], raises ValueError """ - if direction == "top": - T = grid["top_x_first"].values / sizeX * 100 - elif direction == "left": - T = grid["left_y_first"].values / sizeY * 100 - else: - raise ValueError("direction must be either of top or left") - T = T[np.isfinite(T)] - # print(T) - # assert np.all(0 <= T < 100) - - best_model = None - for _ in range(max_stall_count): - init_guess = np.random.uniform(0, 100, size=(3,)) - model = minimize( - compute_inv_liklihood, - init_guess, - args=(T,), - bounds=[(0, 100)] * 3, - method="trust-constr", - ) - if model["x"][0] < prob_uniform_threshold: - if best_model is None: - best_model = model - else: - if np.isclose(best_model["fun"], model["fun"]): - assert len(model["x"]) == 3 - return tuple(map(float, model["x"])) - elif model["fun"] < best_model["fun"]: - best_model = model - assert ( - not best_model is None - ), "Overlap model estimation failed, try raising the value of overlap_prob_uniform_threshold" # noqa : - # best_model_params : Tuple[Float,Float,Float] = tuple(best_model["x"]) - assert len(best_model["x"]) == 3 - return tuple(map(float, best_model["x"])) + translation = np.array( + [ + grid[f"{direction}_y_first"].values / sizeY, + grid[f"{direction}_x_first"].values / sizeX, + ] + ) + translation = translation[:, np.all(np.isfinite(translation), axis=0)] + ifol = IsolationForest() + c = ifol.fit_predict(translation.T) + return tuple(np.median(translation[:, c], axis=1)) def filter_by_overlap_and_correlation( @@ -184,75 +102,13 @@ def filter_outliers(T: pd.Series, isvalid: pd.Series, w: Float = 1.5) -> pd.Seri return isvalid & T.between(q1 - w * iqd, q3 + w * iqd) -def compute_repeatability( - grid: pd.DataFrame, - overlap_n: Float, - overlap_w: Float, - sizeY: Int, - sizeX: Int, - pou: Float, -) -> Tuple[pd.DataFrame, Float]: - """Compute the repeatability of the stage motion. - - Parameters - ---------- - grid : pd.DataFrame - the dataframe for the grid position, with columns "{top|left}_{x|y|ncc}_first" - overlap_n : Float - the estimated overlap for top direction - overlap_w : Float - the estimated overlap for left direction - sizeY : Int - the width of the height images - sizeX : Int - the width of the tile images - pou : Float - the percentile margin for error, by default 3 - - Returns - ------- - grid : pd.DataFrame - the updated dataframe for the grid position - repeatability : Float - the repeatability of the stage motion - """ - grid["top_valid1"] = filter_by_overlap_and_correlation( - grid["top_x_first"], grid["top_ncc_first"], overlap_n, sizeX, pou - ) - grid["left_valid1"] = filter_by_overlap_and_correlation( - grid["left_y_first"], grid["left_ncc_first"], overlap_w, sizeY, pou - ) - grid["top_valid2"] = filter_outliers(grid["top_x_first"], grid["top_valid1"]) - grid["left_valid2"] = filter_outliers(grid["left_y_first"], grid["left_valid1"]) - - if np.any(grid["top_valid2"]): - xs = grid[grid["top_valid2"]]["top_y_first"] - rx_top = np.ceil((xs.max() - xs.min()) / 2) - _, yss = zip(*grid[grid["top_valid2"]].groupby("row")["top_x_first"]) - ry_top = np.ceil(np.max([np.max(ys) - np.min(ys) for ys in yss]) / 2) - r_top = max(rx_top, ry_top) - else: - r_top = 0 # better than failing - - if np.any(grid["left_valid2"]): - ys = grid[grid["left_valid2"]]["left_x_first"] - ry_left = np.ceil((ys.max() - ys.min()) / 2) - _, xss = zip(*grid[grid["left_valid2"]].groupby("col")["left_y_first"]) - rx_left = np.ceil(np.max([np.max(xs) - np.min(xs) for xs in xss]) / 2) - r_left = max(ry_left, rx_left) - else: - r_left = 0 # better than failing - - return grid, max(r_top, r_left) - - def filter_by_repeatability(grid: pd.DataFrame, r: Float) -> pd.DataFrame: """Filter the stage translation by repeatability. Parameters ---------- grid : pd.DataFrame - the dataframe for the grid position, with columns "{top|left}_{x|y|ncc}_first" + the dataframe for the grid position, with columns "{left|top}_{x|y|ncc}_first" r : Float the repeatability value @@ -262,19 +118,7 @@ def filter_by_repeatability(grid: pd.DataFrame, r: Float) -> pd.DataFrame: the updated dataframe for the grid position """ for _, grp in grid.groupby("col"): - isvalid = grp["top_valid2"].astype(bool) - if not any(isvalid): - grid.loc[grp.index, "top_valid3"] = False - else: - medx = grp[isvalid]["top_y_first"].median() - medy = grp[isvalid]["top_x_first"].median() - grid.loc[grp.index, "top_valid3"] = ( - grp["top_y_first"].between(medx - r, medx + r) - & grp["top_x_first"].between(medy - r, medy + r) - & (grp["top_ncc_first"] > 0.5) - ) - for _, grp in grid.groupby("row"): - isvalid = grp["left_valid2"] + isvalid = grp["left_valid2"].astype(bool) if not any(isvalid): grid.loc[grp.index, "left_valid3"] = False else: @@ -285,6 +129,18 @@ def filter_by_repeatability(grid: pd.DataFrame, r: Float) -> pd.DataFrame: & grp["left_x_first"].between(medy - r, medy + r) & (grp["left_ncc_first"] > 0.5) ) + for _, grp in grid.groupby("row"): + isvalid = grp["top_valid2"] + if not any(isvalid): + grid.loc[grp.index, "top_valid3"] = False + else: + medx = grp[isvalid]["top_y_first"].median() + medy = grp[isvalid]["top_x_first"].median() + grid.loc[grp.index, "top_valid3"] = ( + grp["top_y_first"].between(medx - r, medx + r) + & grp["top_x_first"].between(medy - r, medy + r) + & (grp["top_ncc_first"] > 0.5) + ) return grid @@ -295,20 +151,20 @@ def replace_invalid_translations(grid: pd.DataFrame) -> pd.DataFrame: ---------- grid : pd.DataFrame the dataframe for the grid position, - with columns "{top|left}_{x|y}_second" and "{top|left}_valid3" + with columns "{left|top}_{x|y}_second" and "{left|top}_valid3" Returns ------- grid : pd.DataFrame the updatd dataframe for the grid position """ - for direction in ["top", "left"]: + for direction in ["left", "top"]: for key in ["x", "y", "ncc"]: isvalid = grid[f"{direction}_valid3"] grid.loc[isvalid, f"{direction}_{key}_second"] = grid.loc[ isvalid, f"{direction}_{key}_first" ] - for direction, rowcol in zip(["top", "left"], ["col", "row"]): + for direction, rowcol in zip(["left", "top"], ["col", "row"]): for _, grp in grid.groupby(rowcol): isvalid = grp[f"{direction}_valid3"].astype(bool) if any(isvalid): @@ -325,12 +181,12 @@ def replace_invalid_translations(grid: pd.DataFrame) -> pd.DataFrame: f"{direction}_x_first" ].median() grid.loc[grp.index[~isvalid], f"{direction}_ncc_second"] = -1 - for direction, xy in itertools.product(["top", "left"], ["x", "y"]): + for direction, xy in itertools.product(["left", "top"], ["x", "y"]): key = f"{direction}_{xy}_second" isna = pd.isna(grid[key]) grid.loc[isna, key] = grid.loc[~isna, key].median() grid.loc[isna, f"{direction}_ncc_second"] = -1 - for direction, xy in itertools.product(["top", "left"], ["x", "y"]): + for direction, xy in itertools.product(["left", "top"], ["x", "y"]): assert np.all(np.isfinite(grid[f"{direction}_{xy}_second"])) return grid diff --git a/src/m2stitch/_translation_computation.py b/src/m2stitch/_translation_computation.py index 5572314..9fb7a26 100644 --- a/src/m2stitch/_translation_computation.py +++ b/src/m2stitch/_translation_computation.py @@ -164,8 +164,8 @@ def interpret_translation( assert np.all(0 <= xins) and np.all(xins < sizeX) _ncc = -np.infty - x = 0 y = 0 + x = 0 ymagss = [yins, sizeY - yins] ymagss[1][ymagss[0] == 0] = 0 diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index fa8d38a..0ef93c8 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -14,9 +14,10 @@ from ._constrained_refinement import refine_translations from ._global_optimization import compute_final_position from ._global_optimization import compute_maximum_spanning_tree -from ._stage_model import compute_image_overlap -from ._stage_model import compute_repeatability +from ._stage_model import compute_image_overlap2 +from ._stage_model import filter_by_overlap_and_correlation from ._stage_model import filter_by_repeatability +from ._stage_model import filter_outliers from ._stage_model import replace_invalid_translations from ._translation_computation import interpret_translation from ._translation_computation import multi_peak_max @@ -126,17 +127,17 @@ def get_index(col, row): return None grid["top"] = grid.apply( - lambda g: get_index(g["col"] - 1, g["row"]), axis=1 + lambda g: get_index(g["col"], g["row"] - 1), axis=1 ).astype(pd.Int32Dtype()) grid["left"] = grid.apply( - lambda g: get_index(g["col"], g["row"] - 1), axis=1 + lambda g: get_index(g["col"] - 1, g["row"]), axis=1 ).astype(pd.Int32Dtype()) ### dimension order ... m.y.x if position_initial_guess is not None: for j, dimension in enumerate(["y", "x"]): grid[f"{dimension}_pos_init_guess"] = position_initial_guess[:, j] - for direction, dimension in itertools.product(["top", "left"], ["y", "x"]): + for direction, dimension in itertools.product(["left", "top"], ["y", "x"]): for ind, g in grid.iterrows(): i1 = g[direction] if pd.isna(i1): @@ -147,7 +148,7 @@ def get_index(col, row): ) ###### translationComputation ###### - for direction in ["top", "left"]: + for direction in ["left", "top"]: for i2, g in tqdm(grid.iterrows(), total=len(grid)): i1 = g[direction] if pd.isna(i1): @@ -178,25 +179,38 @@ def get_lims(dimension, size): for j, key in enumerate(["ncc", "y", "x"]): grid.loc[i2, f"{direction}_{key}_first"] = max_peak[j] - # top_displacement= - # left_displacement= - prob_uniform_n, mu_n, sigma_n = compute_image_overlap( - grid, "top", sizeY, sizeX, prob_uniform_threshold=overlap_prob_uniform_threshold + left_displacement = compute_image_overlap2( + grid[grid["left_ncc_first"] > 0.5], "left", sizeY, sizeX ) - overlap_n = 100 - mu_n - prob_uniform_w, mu_w, sigma_w = compute_image_overlap( - grid, - "left", - sizeY, - sizeX, - prob_uniform_threshold=overlap_prob_uniform_threshold, + top_displacement = compute_image_overlap2( + grid[grid["top_ncc_first"] > 0.5], "top", sizeY, sizeX ) - overlap_w = 100 - mu_w + overlap_top = np.clip(100 - top_displacement[0] * 100, pou, 100 - pou) + overlap_left = np.clip(100 - left_displacement[1] * 100, pou, 100 - pou) - overlap_n = np.clip(overlap_n, pou, 100 - pou) - overlap_w = np.clip(overlap_w, pou, 100 - pou) + ### compute_repeatability ### + grid["top_valid1"] = filter_by_overlap_and_correlation( + grid["top_y_first"], grid["top_ncc_first"], overlap_top, sizeY, pou + ) + grid["top_valid2"] = filter_outliers(grid["top_y_first"], grid["top_valid1"]) + grid["left_valid1"] = filter_by_overlap_and_correlation( + grid["left_x_first"], grid["left_ncc_first"], overlap_left, sizeX, pou + ) + grid["left_valid2"] = filter_outliers(grid["left_x_first"], grid["left_valid1"]) + + rs = [] + for direction, dims, rowcol in zip(["top", "left"], ["yx", "xy"], ["col", "row"]): + valid_key = f"{direction}_valid2" + valid_grid = grid[grid[valid_key]] + if len(valid_grid) > 0: + w1s = valid_grid[f"{direction}_{dims[0]}_first"] + r1 = np.ceil((w1s.max() - w1s.min()) / 2) + _, w2s = zip(*valid_grid.groupby(rowcol)[f"{direction}_{dims[1]}_first"]) + r2 = np.ceil(np.max([np.max(w2) - np.min(w2) for w2 in w2s]) / 2) + rs.append(max(r1, r2)) + rs.append(0) + r = np.max(rs) - grid, r = compute_repeatability(grid, overlap_n, overlap_w, sizeY, sizeX, pou) grid = filter_by_repeatability(grid, r) grid = replace_invalid_translations(grid) @@ -208,18 +222,8 @@ def get_lims(dimension, size): prop_dict = { "W": sizeY, "H": sizeX, - "overlap_top": overlap_n, - "overlap_left": overlap_w, - "overlap_top_results": { - "prob_uniform": prob_uniform_n, - "mu": mu_n, - "sigma": sigma_n, - }, - "overlap_left_results": { - "prob_uniform": prob_uniform_w, - "mu": mu_w, - "sigma": sigma_w, - }, + "overlap_left": overlap_left, + "overlap_top": overlap_top, "repeatability": r, } if row_col_transpose: diff --git a/tests/test_stitching.py b/tests/test_stitching.py index 5b1f6a4..98c918c 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -25,8 +25,8 @@ def test_stitching(test_image_path: Tuple[npt.NDArray, pd.DataFrame]) -> None: rows = props["row"].to_list() result_df, _ = stitch_images(testimages, rows, cols, row_col_transpose=False) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 2 - 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"])) <= 3 + assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) <= 3 def test_stitching_init_guess( @@ -46,5 +46,5 @@ def test_stitching_init_guess( position_initial_guess=pos_guess, ) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) <= 3 - assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) <= 3 + assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) <= 4 + assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) <= 4 From 9121539629d1818dc3a39049e34929a02f8505af Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:18:35 +0900 Subject: [PATCH 23/31] changed predicter to EllipticEnvelope --- poetry.lock | 154 ++++++++++++++++++++++++++--------- pyproject.toml | 1 + src/m2stitch/_stage_model.py | 10 +-- src/m2stitch/stitching.py | 6 +- 4 files changed, 124 insertions(+), 47 deletions(-) diff --git a/poetry.lock b/poetry.lock index 15b2458..4864d78 100644 --- a/poetry.lock +++ b/poetry.lock @@ -57,7 +57,7 @@ python-versions = ">=3.6.1" [[package]] name = "charset-normalizer" -version = "2.0.11" +version = "2.0.12" description = "The Real First Universal Charset Detector. Open, modern and actively maintained alternative to Chardet." category = "dev" optional = false @@ -172,7 +172,7 @@ python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*" [[package]] name = "importlib-metadata" -version = "4.10.1" +version = "4.11.0" description = "Read metadata from Python packages" category = "dev" optional = false @@ -208,6 +208,14 @@ MarkupSafe = ">=2.0" [package.extras] i18n = ["Babel (>=2.7)"] +[[package]] +name = "joblib" +version = "1.1.0" +description = "Lightweight pipelining with Python functions" +category = "main" +optional = false +python-versions = ">=3.6" + [[package]] name = "livereload" version = "2.6.3" @@ -297,7 +305,7 @@ pyparsing = ">=2.0.2,<3.0.5 || >3.0.5" [[package]] name = "pandas" -version = "1.4.0" +version = "1.4.1" description = "Powerful data structures for data analysis, time series, and statistics" category = "main" optional = false @@ -405,7 +413,7 @@ diagrams = ["jinja2", "railroad-diagrams"] [[package]] name = "pytest" -version = "7.0.0" +version = "7.0.1" description = "pytest: simple powerful testing with Python" category = "dev" optional = false @@ -482,7 +490,7 @@ use_chardet_on_py3 = ["chardet (>=3.0.2,<5)"] [[package]] name = "ruamel.yaml" -version = "0.17.20" +version = "0.17.21" description = "ruamel.yaml is a YAML parser/emitter that supports roundtrip preservation of comments, seq/map flow style, and map key order" category = "dev" optional = false @@ -517,6 +525,26 @@ dparse = ">=0.5.1" packaging = "*" requests = "*" +[[package]] +name = "scikit-learn" +version = "1.0.2" +description = "A set of python modules for machine learning and data mining" +category = "main" +optional = false +python-versions = ">=3.7" + +[package.dependencies] +joblib = ">=0.11" +numpy = ">=1.14.6" +scipy = ">=1.1.0" +threadpoolctl = ">=2.0.0" + +[package.extras] +benchmark = ["matplotlib (>=2.2.3)", "pandas (>=0.25.0)", "memory-profiler (>=0.57.0)"] +docs = ["matplotlib (>=2.2.3)", "scikit-image (>=0.14.5)", "pandas (>=0.25.0)", "seaborn (>=0.9.0)", "memory-profiler (>=0.57.0)", "sphinx (>=4.0.1)", "sphinx-gallery (>=0.7.0)", "numpydoc (>=1.0.0)", "Pillow (>=7.1.2)", "sphinx-prompt (>=1.3.0)", "sphinxext-opengraph (>=0.4.2)"] +examples = ["matplotlib (>=2.2.3)", "scikit-image (>=0.14.5)", "pandas (>=0.25.0)", "seaborn (>=0.9.0)"] +tests = ["matplotlib (>=2.2.3)", "scikit-image (>=0.14.5)", "pandas (>=0.25.0)", "pytest (>=5.0.1)", "pytest-cov (>=2.9.0)", "flake8 (>=3.8.2)", "black (>=21.6b0)", "mypy (>=0.770)", "pyamg (>=4.0.0)"] + [[package]] name = "scipy" version = "1.8.0" @@ -691,6 +719,14 @@ python-versions = ">=3.5" lint = ["flake8", "mypy", "docutils-stubs"] test = ["pytest"] +[[package]] +name = "threadpoolctl" +version = "3.1.0" +description = "threadpoolctl" +category = "main" +optional = false +python-versions = ">=3.6" + [[package]] name = "toml" version = "0.10.2" @@ -745,7 +781,7 @@ test = ["pytest", "typing-extensions", "mypy"] [[package]] name = "typing-extensions" -version = "4.0.1" +version = "4.1.0" description = "Backported and Experimental Type Hints for Python 3.6+" category = "dev" optional = false @@ -817,7 +853,7 @@ testing = ["pytest (>=6)", "pytest-checkdocs (>=2.4)", "pytest-flake8", "pytest- [metadata] lock-version = "1.1" python-versions = ">=3.8,<3.11" -content-hash = "4f771fe50899934788725ecfd1ca152ce39ac4160e0380c4ef7b258001eb687f" +content-hash = "88805c27fb4349cb5f380857ec941ba68659eb7db25e710856e1d79aa120832e" [metadata.files] alabaster = [ @@ -845,8 +881,8 @@ cfgv = [ {file = "cfgv-3.3.1.tar.gz", hash = "sha256:f5a830efb9ce7a445376bb66ec94c638a9787422f96264c98edc6bdeed8ab736"}, ] charset-normalizer = [ - {file = "charset-normalizer-2.0.11.tar.gz", hash = "sha256:98398a9d69ee80548c762ba991a4728bfc3836768ed226b3945908d1a688371c"}, - {file = "charset_normalizer-2.0.11-py3-none-any.whl", hash = "sha256:2842d8f5e82a1f6aa437380934d5e1cd4fcf2003b06fed6940769c164a480a45"}, + {file = "charset-normalizer-2.0.12.tar.gz", hash = "sha256:2857e29ff0d34db842cd7ca3230549d1a697f96ee6d3fb071cfa6c7393832597"}, + {file = "charset_normalizer-2.0.12-py3-none-any.whl", hash = "sha256:6881edbebdb17b39b4eaaa821b438bf6eddffb4468cf344f09f89def34a8b1df"}, ] click = [ {file = "click-8.0.3-py3-none-any.whl", hash = "sha256:353f466495adaeb40b6b5f592f9f91cb22372351c84caeb068132442a4518ef3"}, @@ -928,8 +964,8 @@ imagesize = [ {file = "imagesize-1.3.0.tar.gz", hash = "sha256:cd1750d452385ca327479d45b64d9c7729ecf0b3969a58148298c77092261f9d"}, ] importlib-metadata = [ - {file = "importlib_metadata-4.10.1-py3-none-any.whl", hash = "sha256:899e2a40a8c4a1aec681feef45733de8a6c58f3f6a0dbed2eb6574b4387a77b6"}, - {file = "importlib_metadata-4.10.1.tar.gz", hash = "sha256:951f0d8a5b7260e9db5e41d429285b5f451e928479f19d80818878527d36e95e"}, + {file = "importlib_metadata-4.11.0-py3-none-any.whl", hash = "sha256:6affcdb3aec542dd98df8211e730bba6c5f2bec8288d47bacacde898f548c9ad"}, + {file = "importlib_metadata-4.11.0.tar.gz", hash = "sha256:9e5e553bbba1843cb4a00823014b907616be46ee503d2b9ba001d214a8da218f"}, ] iniconfig = [ {file = "iniconfig-1.1.1-py2.py3-none-any.whl", hash = "sha256:011e24c64b7f47f6ebd835bb12a743f2fbe9a26d4cecaa7f53bc4f35ee9da8b3"}, @@ -939,6 +975,10 @@ jinja2 = [ {file = "Jinja2-3.0.3-py3-none-any.whl", hash = "sha256:077ce6014f7b40d03b47d1f1ca4b0fc8328a692bd284016f806ed0eaca390ad8"}, {file = "Jinja2-3.0.3.tar.gz", hash = "sha256:611bb273cd68f3b993fabdc4064fc858c5b47a973cb5aa7999ec1ba405c87cd7"}, ] +joblib = [ + {file = "joblib-1.1.0-py2.py3-none-any.whl", hash = "sha256:f21f109b3c7ff9d95f8387f752d0d9c34a02aa2f7060c2135f465da0e5160ff6"}, + {file = "joblib-1.1.0.tar.gz", hash = "sha256:4158fcecd13733f8be669be0683b96ebdbbd38d23559f54dca7205aea1bf1e35"}, +] livereload = [ {file = "livereload-2.6.3.tar.gz", hash = "sha256:776f2f865e59fde56490a56bcc6773b6917366bce0c267c60ee8aaf1a0959869"}, ] @@ -1073,27 +1113,27 @@ packaging = [ {file = "packaging-21.3.tar.gz", hash = "sha256:dd47c42927d89ab911e606518907cc2d3a1f38bbd026385970643f9c5b8ecfeb"}, ] pandas = [ - {file = "pandas-1.4.0-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:de62cf699122dcef175988f0714678e59c453dc234c5b47b7136bfd7641e3c8c"}, - {file = "pandas-1.4.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:46a18572f3e1cb75db59d9461940e9ba7ee38967fa48dd58f4139197f6e32280"}, - {file = "pandas-1.4.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:73f7da2ccc38cc988b74e5400b430b7905db5f2c413ff215506bea034eaf832d"}, - {file = "pandas-1.4.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:5229c95db3a907451dacebc551492db6f7d01743e49bbc862f4a6010c227d187"}, - {file = "pandas-1.4.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:fe454180ad31bbbe1e5d111b44443258730467f035e26b4e354655ab59405871"}, - {file = "pandas-1.4.0-cp310-cp310-win_amd64.whl", hash = "sha256:784cca3f69cfd7f6bd7c7fdb44f2bbab17e6de55725e9ff36d6f382510dfefb5"}, - {file = "pandas-1.4.0-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:de8f8999864399529e8514a2e6bfe00fd161f0a667903655552ed12e583ae3cb"}, - {file = "pandas-1.4.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:0f19504f2783526fb5b4de675ea69d68974e21c1624f4b92295d057a31d5ec5f"}, - {file = "pandas-1.4.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:f045bb5c6bfaba536089573bf97d6b8ccc7159d951fe63904c395a5e486fbe14"}, - {file = "pandas-1.4.0-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:5280d057ddae06fe4a3cd6aa79040b8c205cd6dd21743004cf8635f39ed01712"}, - {file = "pandas-1.4.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:1f3b74335390dda49f5d5089fab71958812bf56f42aa27663ee4c16d19f4f1c5"}, - {file = "pandas-1.4.0-cp38-cp38-win32.whl", hash = "sha256:51e5da3802aaee1aa4254108ffaf1129a15fb3810b7ce8da1ec217c655b418f5"}, - {file = "pandas-1.4.0-cp38-cp38-win_amd64.whl", hash = "sha256:f103a5cdcd66cb18882ccdc18a130c31c3cfe3529732e7f10a8ab3559164819c"}, - {file = "pandas-1.4.0-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:4a8d5a200f8685e7ea562b2f022c77ab7cb82c1ca5b240e6965faa6f84e5c1e9"}, - {file = "pandas-1.4.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:b5af258c7b090cca7b742cf2bd67ad1919aa9e4e681007366c9edad2d6a3d42b"}, - {file = "pandas-1.4.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:156aac90dd7b303bf0b91bae96c0503212777f86c731e41929c571125d26c8e9"}, - {file = "pandas-1.4.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:2dad075089e17a72391de33021ad93720aff258c3c4b68c78e1cafce7e447045"}, - {file = "pandas-1.4.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:1d59c958d6b8f96fdf850c7821571782168d5acfe75ccf78cd8d1ac15fb921df"}, - {file = "pandas-1.4.0-cp39-cp39-win32.whl", hash = "sha256:55ec0e192eefa26d823fc25a1f213d6c304a3592915f368e360652994cdb8d9a"}, - {file = "pandas-1.4.0-cp39-cp39-win_amd64.whl", hash = "sha256:23c04dab11f3c6359cfa7afa83d3d054a8f8c283d773451184d98119ef54da97"}, - {file = "pandas-1.4.0.tar.gz", hash = "sha256:cdd76254c7f0a1583bd4e4781fb450d0ebf392e10d3f12e92c95575942e37df5"}, + {file = "pandas-1.4.1-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:3dfb32ed50122fe8c5e7f2b8d97387edd742cc78f9ec36f007ee126cd3720907"}, + {file = "pandas-1.4.1-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:0259cd11e7e6125aaea3af823b80444f3adad6149ff4c97fef760093598b3e34"}, + {file = "pandas-1.4.1-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:96e9ece5759f9b47ae43794b6359bbc54805d76e573b161ae770c1ea59393106"}, + {file = "pandas-1.4.1-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:508c99debccd15790d526ce6b1624b97a5e1e4ca5b871319fb0ebfd46b8f4dad"}, + {file = "pandas-1.4.1-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:e6a7bbbb7950063bfc942f8794bc3e31697c020a14f1cd8905fc1d28ec674a01"}, + {file = "pandas-1.4.1-cp310-cp310-win_amd64.whl", hash = "sha256:c614001129b2a5add5e3677c3a213a9e6fd376204cb8d17c04e84ff7dfc02a73"}, + {file = "pandas-1.4.1-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:4e1176f45981c8ccc8161bc036916c004ca51037a7ed73f2d2a9857e6dbe654f"}, + {file = "pandas-1.4.1-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:bbb15ad79050e8b8d39ec40dd96a30cd09b886a2ae8848d0df1abba4d5502a67"}, + {file = "pandas-1.4.1-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:6d6ad1da00c7cc7d8dd1559a6ba59ba3973be6b15722d49738b2be0977eb8a0c"}, + {file = "pandas-1.4.1-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:358b0bc98a5ff067132d23bf7a2242ee95db9ea5b7bbc401cf79205f11502fd3"}, + {file = "pandas-1.4.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:6105af6533f8b63a43ea9f08a2ede04e8f43e49daef0209ab0d30352bcf08bee"}, + {file = "pandas-1.4.1-cp38-cp38-win32.whl", hash = "sha256:04dd15d9db538470900c851498e532ef28d4e56bfe72c9523acb32042de43dfb"}, + {file = "pandas-1.4.1-cp38-cp38-win_amd64.whl", hash = "sha256:1b384516dbb4e6aae30e3464c2e77c563da5980440fbdfbd0968e3942f8f9d70"}, + {file = "pandas-1.4.1-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:f02e85e6d832be37d7f16cf6ac8bb26b519ace3e5f3235564a91c7f658ab2a43"}, + {file = "pandas-1.4.1-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:0b1a13f647e4209ed7dbb5da3497891d0045da9785327530ab696417ef478f84"}, + {file = "pandas-1.4.1-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:19f7c632436b1b4f84615c3b127bbd7bc603db95e3d4332ed259dc815c9aaa26"}, + {file = "pandas-1.4.1-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:7ea47ba1d6f359680130bd29af497333be6110de8f4c35b9211eec5a5a9630fa"}, + {file = "pandas-1.4.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:2e5a7a1e0ecaac652326af627a3eca84886da9e667d68286866d4e33f6547caf"}, + {file = "pandas-1.4.1-cp39-cp39-win32.whl", hash = "sha256:1d85d5f6be66dfd6d1d8d13b9535e342a2214260f1852654b19fa4d7b8d1218b"}, + {file = "pandas-1.4.1-cp39-cp39-win_amd64.whl", hash = "sha256:3129a35d9dad1d80c234dd78f8f03141b914395d23f97cf92a366dcd19f8f8bf"}, + {file = "pandas-1.4.1.tar.gz", hash = "sha256:8db93ec98ac7cb5f8ac1420c10f5e3c43533153f253fe7fb6d891cf5aa2b80d2"}, ] pandas-stubs = [ {file = "pandas-stubs-1.2.0.47.tar.gz", hash = "sha256:78801efb504abe75c2cf73cb3469b97003080e80fb03aa2a406ca4928f3fd887"}, @@ -1128,8 +1168,8 @@ pyparsing = [ {file = "pyparsing-3.0.7.tar.gz", hash = "sha256:18ee9022775d270c55187733956460083db60b37d0d0fb357445f3094eed3eea"}, ] pytest = [ - {file = "pytest-7.0.0-py3-none-any.whl", hash = "sha256:42901e6bd4bd4a0e533358a86e848427a49005a3256f657c5c8f8dd35ef137a9"}, - {file = "pytest-7.0.0.tar.gz", hash = "sha256:dad48ffda394e5ad9aa3b7d7ddf339ed502e5e365b1350e0af65f4a602344b11"}, + {file = "pytest-7.0.1-py3-none-any.whl", hash = "sha256:9ce3ff477af913ecf6321fe337b93a2c0dcf2a0a1439c43f5452112c1e4280db"}, + {file = "pytest-7.0.1.tar.gz", hash = "sha256:e30905a0c131d3d94b89624a1cc5afec3e0ba2fbdb151867d8e0ebd49850f171"}, ] pytest-datadir = [ {file = "pytest-datadir-1.3.1.tar.gz", hash = "sha256:d3af1e738df87515ee509d6135780f25a15959766d9c2b2dbe02bf4fb979cb18"}, @@ -1183,8 +1223,8 @@ requests = [ {file = "requests-2.27.1.tar.gz", hash = "sha256:68d7c56fd5a8999887728ef304a6d12edc7be74f1cfa47714fc8b414525c9a61"}, ] "ruamel.yaml" = [ - {file = "ruamel.yaml-0.17.20-py3-none-any.whl", hash = "sha256:810eef9c46523a3f77479c66267a4708255ebe806a2d540078408c2227f011af"}, - {file = "ruamel.yaml-0.17.20.tar.gz", hash = "sha256:4b8a33c1efb2b443a93fcaafcfa4d2e445f8e8c29c528d9f5cdafb7cc9e4004c"}, + {file = "ruamel.yaml-0.17.21-py3-none-any.whl", hash = "sha256:742b35d3d665023981bd6d16b3d24248ce5df75fdb4e2924e93a05c1f8b61ca7"}, + {file = "ruamel.yaml-0.17.21.tar.gz", hash = "sha256:8b7ce697a2f212752a35c1ac414471dc16c424c9573be4926b56ff3f5d23b7af"}, ] "ruamel.yaml.clib" = [ {file = "ruamel.yaml.clib-0.2.6-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:6e7be2c5bcb297f5b82fee9c665eb2eb7001d1050deaba8471842979293a80b0"}, @@ -1217,6 +1257,40 @@ safety = [ {file = "safety-1.10.3-py2.py3-none-any.whl", hash = "sha256:5f802ad5df5614f9622d8d71fedec2757099705c2356f862847c58c6dfe13e84"}, {file = "safety-1.10.3.tar.gz", hash = "sha256:30e394d02a20ac49b7f65292d19d38fa927a8f9582cdfd3ad1adbbc66c641ad5"}, ] +scikit-learn = [ + {file = "scikit-learn-1.0.2.tar.gz", hash = "sha256:b5870959a5484b614f26d31ca4c17524b1b0317522199dc985c3b4256e030767"}, + {file = "scikit_learn-1.0.2-cp310-cp310-macosx_10_13_x86_64.whl", hash = "sha256:da3c84694ff693b5b3194d8752ccf935a665b8b5edc33a283122f4273ca3e687"}, + {file = "scikit_learn-1.0.2-cp310-cp310-macosx_12_0_arm64.whl", hash = "sha256:75307d9ea39236cad7eea87143155eea24d48f93f3a2f9389c817f7019f00705"}, + {file = "scikit_learn-1.0.2-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:f14517e174bd7332f1cca2c959e704696a5e0ba246eb8763e6c24876d8710049"}, + {file = "scikit_learn-1.0.2-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d9aac97e57c196206179f674f09bc6bffcd0284e2ba95b7fe0b402ac3f986023"}, + {file = "scikit_learn-1.0.2-cp310-cp310-win_amd64.whl", hash = "sha256:d93d4c28370aea8a7cbf6015e8a669cd5d69f856cc2aa44e7a590fb805bb5583"}, + {file = "scikit_learn-1.0.2-cp37-cp37m-macosx_10_13_x86_64.whl", hash = "sha256:85260fb430b795d806251dd3bb05e6f48cdc777ac31f2bcf2bc8bbed3270a8f5"}, + {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:a053a6a527c87c5c4fa7bf1ab2556fa16d8345cf99b6c5a19030a4a7cd8fd2c0"}, + {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:245c9b5a67445f6f044411e16a93a554edc1efdcce94d3fc0bc6a4b9ac30b752"}, + {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:158faf30684c92a78e12da19c73feff9641a928a8024b4fa5ec11d583f3d8a87"}, + {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:08ef968f6b72033c16c479c966bf37ccd49b06ea91b765e1cc27afefe723920b"}, + {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:16455ace947d8d9e5391435c2977178d0ff03a261571e67f627c8fee0f9d431a"}, + {file = "scikit_learn-1.0.2-cp37-cp37m-win32.whl", hash = "sha256:2f3b453e0b149898577e301d27e098dfe1a36943f7bb0ad704d1e548efc3b448"}, + {file = "scikit_learn-1.0.2-cp37-cp37m-win_amd64.whl", hash = "sha256:46f431ec59dead665e1370314dbebc99ead05e1c0a9df42f22d6a0e00044820f"}, + {file = "scikit_learn-1.0.2-cp38-cp38-macosx_10_13_x86_64.whl", hash = "sha256:ff3fa8ea0e09e38677762afc6e14cad77b5e125b0ea70c9bba1992f02c93b028"}, + {file = "scikit_learn-1.0.2-cp38-cp38-macosx_12_0_arm64.whl", hash = "sha256:9369b030e155f8188743eb4893ac17a27f81d28a884af460870c7c072f114243"}, + {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:7d6b2475f1c23a698b48515217eb26b45a6598c7b1840ba23b3c5acece658dbb"}, + {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:285db0352e635b9e3392b0b426bc48c3b485512d3b4ac3c7a44ec2a2ba061e66"}, + {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:5cb33fe1dc6f73dc19e67b264dbb5dde2a0539b986435fdd78ed978c14654830"}, + {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:b1391d1a6e2268485a63c3073111fe3ba6ec5145fc957481cfd0652be571226d"}, + {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:bc3744dabc56b50bec73624aeca02e0def06b03cb287de26836e730659c5d29c"}, + {file = "scikit_learn-1.0.2-cp38-cp38-win32.whl", hash = "sha256:a999c9f02ff9570c783069f1074f06fe7386ec65b84c983db5aeb8144356a355"}, + {file = "scikit_learn-1.0.2-cp38-cp38-win_amd64.whl", hash = "sha256:7626a34eabbf370a638f32d1a3ad50526844ba58d63e3ab81ba91e2a7c6d037e"}, + {file = "scikit_learn-1.0.2-cp39-cp39-macosx_10_13_x86_64.whl", hash = "sha256:a90b60048f9ffdd962d2ad2fb16367a87ac34d76e02550968719eb7b5716fd10"}, + {file = "scikit_learn-1.0.2-cp39-cp39-macosx_12_0_arm64.whl", hash = "sha256:7a93c1292799620df90348800d5ac06f3794c1316ca247525fa31169f6d25855"}, + {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:eabceab574f471de0b0eb3f2ecf2eee9f10b3106570481d007ed1c84ebf6d6a1"}, + {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:55f2f3a8414e14fbee03782f9fe16cca0f141d639d2b1c1a36779fa069e1db57"}, + {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:80095a1e4b93bd33261ef03b9bc86d6db649f988ea4dbcf7110d0cded8d7213d"}, + {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:fa38a1b9b38ae1fad2863eff5e0d69608567453fdfc850c992e6e47eb764e846"}, + {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ff746a69ff2ef25f62b36338c615dd15954ddc3ab8e73530237dd73235e76d62"}, + {file = "scikit_learn-1.0.2-cp39-cp39-win32.whl", hash = "sha256:e174242caecb11e4abf169342641778f68e1bfaba80cd18acd6bc84286b9a534"}, + {file = "scikit_learn-1.0.2-cp39-cp39-win_amd64.whl", hash = "sha256:b54a62c6e318ddbfa7d22c383466d38d2ee770ebdb5ddb668d56a099f6eaf75f"}, +] scipy = [ {file = "scipy-1.8.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:87b01c7d5761e8a266a0fbdb9d88dcba0910d63c1c671bdb4d99d29f469e9e03"}, {file = "scipy-1.8.0-cp310-cp310-macosx_12_0_arm64.whl", hash = "sha256:ae3e327da323d82e918e593460e23babdce40d7ab21490ddf9fc06dec6b91a18"}, @@ -1290,6 +1364,10 @@ sphinxcontrib-serializinghtml = [ {file = "sphinxcontrib-serializinghtml-1.1.5.tar.gz", hash = "sha256:aa5f6de5dfdf809ef505c4895e51ef5c9eac17d0f287933eb49ec495280b6952"}, {file = "sphinxcontrib_serializinghtml-1.1.5-py2.py3-none-any.whl", hash = "sha256:352a9a00ae864471d3a7ead8d7d79f5fc0b57e8b3f95e9867eb9eb28999b92fd"}, ] +threadpoolctl = [ + {file = "threadpoolctl-3.1.0-py3-none-any.whl", hash = "sha256:8b99adda265feb6773280df41eece7b2e6561b772d21ffd52e372f999024907b"}, + {file = "threadpoolctl-3.1.0.tar.gz", hash = "sha256:a335baacfaa4400ae1f0d8e3a58d6674d2f8828e3716bb2802c44955ad391380"}, +] toml = [ {file = "toml-0.10.2-py2.py3-none-any.whl", hash = "sha256:806143ae5bfb6a3c6e736a764057db0e6a0e05e338b5630894a5f779cabb4f9b"}, {file = "toml-0.10.2.tar.gz", hash = "sha256:b3bda1d108d5dd99f4a20d24d9c348e91c4db7ab1b749200bded2f839ccbe68f"}, @@ -1350,8 +1428,8 @@ typeguard = [ {file = "typeguard-2.13.3.tar.gz", hash = "sha256:00edaa8da3a133674796cf5ea87d9f4b4c367d77476e185e80251cc13dfbb8c4"}, ] typing-extensions = [ - {file = "typing_extensions-4.0.1-py3-none-any.whl", hash = "sha256:7f001e5ac290a0c0401508864c7ec868be4e701886d5b573a9528ed3973d9d3b"}, - {file = "typing_extensions-4.0.1.tar.gz", hash = "sha256:4ca091dea149f945ec56afb48dae714f21e8692ef22a395223bcd328961b6a0e"}, + {file = "typing_extensions-4.1.0-py3-none-any.whl", hash = "sha256:c13180fbaa7cd97065a4915ceba012bdb31dc34743e63ddee16360161d358414"}, + {file = "typing_extensions-4.1.0.tar.gz", hash = "sha256:ba97c5143e5bb067b57793c726dd857b1671d4b02ced273ca0538e71ff009095"}, ] urllib3 = [ {file = "urllib3-1.26.8-py2.py3-none-any.whl", hash = "sha256:000ca7f471a233c2251c6c7023ee85305721bfdf18621ebff4fd17a8653427ed"}, diff --git a/pyproject.toml b/pyproject.toml index 4bd9880..108a08e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,6 +32,7 @@ pandas-stubs = "^1.2.0" networkx = "^2.5.1" scipy = "^1.6.3" tqdm = "^4.60.0" +scikit-learn = "^1.0.2" [tool.poetry.dev-dependencies] pytest = "^7.0.0" diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index a48784c..f8c832f 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -1,19 +1,16 @@ import itertools +from typing import Callable from typing import Tuple import numpy as np import pandas as pd -from sklearn.ensemble import IsolationForest from ._typing_utils import Float from ._typing_utils import Int def compute_image_overlap2( - grid: pd.DataFrame, - direction: str, - sizeY: Int, - sizeX: Int, + grid: pd.DataFrame, direction: str, sizeY: Int, sizeX: Int, predictor: Callable ) -> Tuple[Float, Float]: """Compute the value of the image overlap. @@ -45,8 +42,7 @@ def compute_image_overlap2( ] ) translation = translation[:, np.all(np.isfinite(translation), axis=0)] - ifol = IsolationForest() - c = ifol.fit_predict(translation.T) + c = predictor.fit_predict(translation.T) return tuple(np.median(translation[:, c], axis=1)) diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index 0ef93c8..6c2191a 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -9,6 +9,7 @@ import numpy as np import pandas as pd +from sklearn.covariance import EllipticEnvelope from tqdm import tqdm from ._constrained_refinement import refine_translations @@ -179,11 +180,12 @@ def get_lims(dimension, size): for j, key in enumerate(["ncc", "y", "x"]): grid.loc[i2, f"{direction}_{key}_first"] = max_peak[j] + predictor = EllipticEnvelope(contamination=0.4) left_displacement = compute_image_overlap2( - grid[grid["left_ncc_first"] > 0.5], "left", sizeY, sizeX + grid[grid["left_ncc_first"] > 0.5], "left", sizeY, sizeX, predictor ) top_displacement = compute_image_overlap2( - grid[grid["top_ncc_first"] > 0.5], "top", sizeY, sizeX + grid[grid["top_ncc_first"] > 0.5], "top", sizeY, sizeX, predictor ) overlap_top = np.clip(100 - top_displacement[0] * 100, pou, 100 - pou) overlap_left = np.clip(100 - left_displacement[1] * 100, pou, 100 - pou) From 81d9171b6d5ff293d5ac319dc047c70a8af0bdb5 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:25:06 +0900 Subject: [PATCH 24/31] try to solve mypy problem --- src/m2stitch/_stage_model.py | 2 +- tests/data/testimages_props.csv | 2 +- tests/test_stitching.py | 8 ++++---- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index f8c832f..188cd71 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -43,7 +43,7 @@ def compute_image_overlap2( ) translation = translation[:, np.all(np.isfinite(translation), axis=0)] c = predictor.fit_predict(translation.T) - return tuple(np.median(translation[:, c], axis=1)) + return tuple(np.median(translation[:, c], axis=1).astype(np.float32)) def filter_by_overlap_and_correlation( diff --git a/tests/data/testimages_props.csv b/tests/data/testimages_props.csv index fb30292..16bd132 100644 --- a/tests/data/testimages_props.csv +++ b/tests/data/testimages_props.csv @@ -3,7 +3,7 @@ 1,3,0,0,5000 2,3,1,1349,5001 3,2,1,1350,3331 -4,1,1,1353,1664 +4,1,1,1353,1668 5,0,2,2703,0 6,1,2,2701,1664 7,2,2,2700,3331 diff --git a/tests/test_stitching.py b/tests/test_stitching.py index 98c918c..8a7be07 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -25,8 +25,8 @@ def test_stitching(test_image_path: Tuple[npt.NDArray, pd.DataFrame]) -> None: rows = props["row"].to_list() result_df, _ = stitch_images(testimages, rows, cols, row_col_transpose=False) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) <= 3 - assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) <= 3 + assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 2 + assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 2 def test_stitching_init_guess( @@ -46,5 +46,5 @@ def test_stitching_init_guess( position_initial_guess=pos_guess, ) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) <= 4 - assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) <= 4 + assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 2 + assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 2 From b1e1d4f773efb559695e639aff8d84df60ae2ffa Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:29:32 +0900 Subject: [PATCH 25/31] solving mypy ... --- src/m2stitch/_stage_model.py | 4 ++-- tests/test_stitching.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index 188cd71..a08be68 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -1,16 +1,16 @@ import itertools -from typing import Callable from typing import Tuple import numpy as np import pandas as pd +from sklearn.base import BaseEstimator from ._typing_utils import Float from ._typing_utils import Int def compute_image_overlap2( - grid: pd.DataFrame, direction: str, sizeY: Int, sizeX: Int, predictor: Callable + grid: pd.DataFrame, direction: str, sizeY: Int, sizeX: Int, predictor: BaseEstimator ) -> Tuple[Float, Float]: """Compute the value of the image overlap. diff --git a/tests/test_stitching.py b/tests/test_stitching.py index 8a7be07..70de5d0 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -25,8 +25,8 @@ def test_stitching(test_image_path: Tuple[npt.NDArray, pd.DataFrame]) -> None: rows = props["row"].to_list() result_df, _ = stitch_images(testimages, rows, cols, row_col_transpose=False) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 2 - 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"])) < 3 + assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 3 def test_stitching_init_guess( @@ -46,5 +46,5 @@ def test_stitching_init_guess( position_initial_guess=pos_guess, ) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 2 - 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"])) < 3 + assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 3 From b09187a3a65da6dfe08a6f34eeda90a5b0d54d5e Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:34:58 +0900 Subject: [PATCH 26/31] changed threshold for tests --- tests/data/testimages_props.csv | 32 ++++++++++++++++---------------- tests/test_stitching.py | 8 ++++---- 2 files changed, 20 insertions(+), 20 deletions(-) diff --git a/tests/data/testimages_props.csv b/tests/data/testimages_props.csv index 16bd132..4221a0c 100644 --- a/tests/data/testimages_props.csv +++ b/tests/data/testimages_props.csv @@ -1,16 +1,16 @@ -,col,row,y_pos,x_pos -0,2,0,1,3330 -1,3,0,0,5000 -2,3,1,1349,5001 -3,2,1,1350,3331 -4,1,1,1353,1668 -5,0,2,2703,0 -6,1,2,2701,1664 -7,2,2,2700,3331 -8,3,2,2699,5002 -9,3,3,4049,5003 -10,2,3,4050,3333 -11,1,3,4051,1666 -12,0,3,4052,1 -13,0,4,5401,1 -14,1,4,5400,1666 +,col,row,y_pos,x_pos,allowed_error +0,2,0,1,3330,1 +1,3,0,0,5000,1 +2,3,1,1349,5001,1 +3,2,1,1350,3331,1 +4,1,1,1353,1668,4 +5,0,2,2703,0,1 +6,1,2,2701,1664,1 +7,2,2,2700,3331,1 +8,3,2,2699,5002,1 +9,3,3,4049,5003,1 +10,2,3,4050,3333,1 +11,1,3,4051,1666,1 +12,0,3,4052,1,1 +13,0,4,5401,1,1 +14,1,4,5400,1666,1 diff --git a/tests/test_stitching.py b/tests/test_stitching.py index 70de5d0..8de11c4 100644 --- a/tests/test_stitching.py +++ b/tests/test_stitching.py @@ -25,8 +25,8 @@ def test_stitching(test_image_path: Tuple[npt.NDArray, pd.DataFrame]) -> None: rows = props["row"].to_list() result_df, _ = stitch_images(testimages, rows, cols, row_col_transpose=False) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 3 - assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 3 + assert all(np.abs(result_df["y_pos"] - props["y_pos"]) <= props["allowed_error"]) + assert all(np.abs(result_df["x_pos"] - props["x_pos"]) <= props["allowed_error"]) def test_stitching_init_guess( @@ -46,5 +46,5 @@ def test_stitching_init_guess( position_initial_guess=pos_guess, ) assert np.array_equal(result_df.index, props.index) - assert np.max(np.abs(result_df["y_pos"] - props["y_pos"])) < 3 - assert np.max(np.abs(result_df["x_pos"] - props["x_pos"])) < 3 + assert all(np.abs(result_df["y_pos"] - props["y_pos"]) <= props["allowed_error"]) + assert all(np.abs(result_df["x_pos"] - props["x_pos"]) <= props["allowed_error"]) From ece4d33c3a87d95e285643b31997a6477de4a11a Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:39:22 +0900 Subject: [PATCH 27/31] mypy problem solved --- src/m2stitch/_stage_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index a08be68..b53a3ad 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -43,7 +43,7 @@ def compute_image_overlap2( ) translation = translation[:, np.all(np.isfinite(translation), axis=0)] c = predictor.fit_predict(translation.T) - return tuple(np.median(translation[:, c], axis=1).astype(np.float32)) + return tuple(map(float, np.median(translation[:, c], axis=1))) def filter_by_overlap_and_correlation( From 802d153c70269690d9335fed514bc1e7167b18b5 Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:40:07 +0900 Subject: [PATCH 28/31] bumped version to 0.4.0 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 108a08e..b54c189 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "m2stitch" -version = "0.3.0" +version = "0.4.0" description = "M2Stitch" authors = ["Yohsuke Fukai "] license = "MIT" From 88105a26fad61b08b9f8c464e4bba559697d551a Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:41:38 +0900 Subject: [PATCH 29/31] added comment --- src/m2stitch/stitching.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/m2stitch/stitching.py b/src/m2stitch/stitching.py index 6c2191a..26bcaae 100644 --- a/src/m2stitch/stitching.py +++ b/src/m2stitch/stitching.py @@ -181,6 +181,7 @@ def get_lims(dimension, size): grid.loc[i2, f"{direction}_{key}_first"] = max_peak[j] predictor = EllipticEnvelope(contamination=0.4) + # TODO make threshold adjustable:w left_displacement = compute_image_overlap2( grid[grid["left_ncc_first"] > 0.5], "left", sizeY, sizeX, predictor ) From 6cbdfc4261caffbb8f9b1fa4fc28eea01f4771ed Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:44:37 +0900 Subject: [PATCH 30/31] trying to solve mypy things... --- src/m2stitch/_stage_model.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index b53a3ad..e6b314a 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -43,7 +43,9 @@ def compute_image_overlap2( ) translation = translation[:, np.all(np.isfinite(translation), axis=0)] c = predictor.fit_predict(translation.T) - return tuple(map(float, np.median(translation[:, c], axis=1))) + res = np.median(translation[:, c], axis=1) + assert len(res) == 2 + return tuple(res) def filter_by_overlap_and_correlation( From 444c5173174b5f5dff55e8892cd9e10bf210bacb Mon Sep 17 00:00:00 2001 From: Yohsuke Fukai Date: Sun, 13 Feb 2022 16:46:08 +0900 Subject: [PATCH 31/31] hopefully solved mypy --- src/m2stitch/_stage_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/m2stitch/_stage_model.py b/src/m2stitch/_stage_model.py index e6b314a..ec409f0 100644 --- a/src/m2stitch/_stage_model.py +++ b/src/m2stitch/_stage_model.py @@ -11,7 +11,7 @@ def compute_image_overlap2( grid: pd.DataFrame, direction: str, sizeY: Int, sizeX: Int, predictor: BaseEstimator -) -> Tuple[Float, Float]: +) -> Tuple[Float, ...]: """Compute the value of the image overlap. Parameters