Commit e2a6ee31 authored by Sean Fitzgibbon's avatar Sean Fitzgibbon
Browse files

Initial commit of new optimised_xfm algorithm

parent b2314e38
......@@ -48,6 +48,19 @@ def register_chart_to_slide(
do_plots: Optional[bool] = None,
justify: Optional[str] = None,
):
'''
It takes a chart and a slide and aligns the chart to the slide.
Args:
chart (str): The path to the chart file.
slide (str): The slide image.
slide_res (float): The resolution of the slide.
outdir (str): The output directory for the results.
boundary_key (Optional[str]): The name of the boundary contour in the chart.
config (Optional[str]): Path to the YAML configuration file.
do_plots (Optional[bool]): Whether to plot the results of the alignment.
justify (Optional[str]): str, optional
'''
if config is None:
config = util.get_resource("chart.yaml")
......@@ -110,16 +123,19 @@ def register_chart_to_slide(
np.savetxt(f"{outdir}/chart-to-image-init.xfm", init_xfm.params)
# refine alignment (to mask edges)
opt_xfm = optimise_xfm(img, slide_res, edge_crds_cat, init_xfm)
# refine alignment (to mask edges)
# print(opt_xfm)
mask = img_props['mask']
opt_xfm, mdist = optimise_xfm(mask, slide_res, edge_crds_cat[::3,:], init_xfm)
tr_x, tr_y = opt_xfm.translation
print(
f"Optimised XFM - Rotation: {opt_xfm.rotation:0.5f}, Translation: [{tr_x:0.5f} {tr_y:0.5f}], Scale: {opt_xfm.scale:0.5f}"
)
plt.plot(mdist)
plt.show()
np.savetxt(f"{outdir}/chart-to-image.xfm", opt_xfm.params)
# apply optimised xfm to contours and cells and save
......@@ -152,6 +168,9 @@ def register_chart_to_slide(
if do_plots:
# TODO: plot normals
# TODO: plot iteration errors for optimisation
fig, ax = plt.subplots(2, 2, figsize=(40, 40))
ax = np.ravel(ax)
......@@ -316,7 +335,22 @@ def segment_foreground(
return brainmask
def optimise_xfm(image, image_resolution, contour, xfm_init):
def optimise_xfm(mask, mask_resolution, contour, xfm, n_iter=100):
mdist = np.zeros(n_iter)
# TODO: move line extent values to config
max_line_extent = 1.5
min_line_extent = 0.01
for idx in range(n_iter):
xfm, mdist[idx] = optimise_xfm_worker(mask, mask_resolution, contour, xfm, line_extent=max_line_extent)
max_line_extent = np.maximum(max_line_extent*0.9, min_line_extent)
return xfm, mdist
def optimise_xfm_worker(image, image_resolution, contour, xfm_init, line_extent=1.5):
"""
Refine contour by sampling image along normal (to each edge point) and looking for big step change.
Calculate a new transform (xfm) using refined contour points.
......@@ -332,41 +366,61 @@ def optimise_xfm(image, image_resolution, contour, xfm_init):
# TODO: move these line extents to the config file
# line_smpl = np.linspace(-0.03, 0.15, 20)
line_smpl = np.linspace(-0.2, 0.2, 20)
line_smpl = np.linspace(-line_extent, line_extent, 1000)
line_x = edge_init_x[:, np.newaxis] + nrml_x[:, np.newaxis] * line_smpl
line_y = edge_init_y[:, np.newaxis] + nrml_y[:, np.newaxis] * line_smpl
brainmask = segment_foreground(image)
# brainmask = segment_foreground(image)
brainmask = image
# if do_plots:
# fig, ax = plt.subplots(figsize=(20, 30))
# sample image along normal line
line_x = np.round(line_x / image_resolution).astype(int)
line_y = np.round(line_y / image_resolution).astype(int)
# extent = np.array([0, img.shape[1], img.shape[0], 0]) * img_res
# ax.imshow(brainmask, extent=extent, cmap='gray')
y, x = brainmask.shape
line_int = brainmask[np.clip(line_y, 0, y - 1), np.clip(line_x, 0, x - 1)]
# for line_x0, line_y0 in zip(line_x, line_y):
# ax.plot(line_x0, line_y0, "b.")
# exclude normals that cross the mask edge twice
single_crossing_idx = np.sum(np.diff(line_int, axis=1), axis=1)==1
debug = False
if debug:
np.savez(
f"./debug.npz",
line_x=line_x,
line_y=line_y,
line_int=line_int,
brainmask=brainmask,
image=image,
image_resolution=image_resolution,
contour=contour,
normals=normals,
single_crossing_idx=single_crossing_idx,
edge_coords_init=edge_coords_init,
)
fig, ax = plt.subplots(figsize=(20, 30))
# ax.plot(edge_init_x, edge_init_y, "r.")
extent = np.array([0, image.shape[1], image.shape[0], 0])
ax.imshow(brainmask, extent=extent, cmap="gray")
# ax.set_xlabel("mm")
# ax.set_ylabel("mm")
# ax.set_title("Brainmask + normals")
for line_x0, line_y0 in zip(line_x[single_crossing_idx, :], line_y[single_crossing_idx, :]):
ax.plot(line_x0, line_y0, "b.")
# plt.show()
# # fig.savefig(f"{OUTDIR}/normals.png", bbox_inches="tight", dpi=300)
ax.plot(edge_init_x, edge_init_y, "r.")
# sample image along normal line
ax.set_xlabel("mm")
ax.set_ylabel("mm")
ax.set_title("Brainmask + normals")
line_x = np.round(line_x / image_resolution).astype(int)
line_y = np.round(line_y / image_resolution).astype(int)
plt.show()
y, x = brainmask.shape
line_int = brainmask[np.clip(line_y, 0, y - 1), np.clip(line_x, 0, x - 1)]
# exlude constant rows
constant_idx = np.all(line_int == line_int[:, 0][:, np.newaxis], axis=1)
# constant_idx = np.all(line_int == line_int[:, 0][:, np.newaxis], axis=1)
min_idx = np.argmax(np.abs(np.diff(line_int, axis=-1)), axis=-1)
......@@ -379,11 +433,16 @@ def optimise_xfm(image, image_resolution, contour, xfm_init):
)
opt_xfm = transform.SimilarityTransform()
opt_xfm.estimate(
contour[~constant_idx, :], refined_edge_coords[~constant_idx, :]
)
opt_xfm.estimate(contour[single_crossing_idx, :], refined_edge_coords[single_crossing_idx, :])
# calculate error (mean distance)
opt_edge_coords = apply_xfm(opt_xfm, contour)
dist = np.sqrt(np.sum((opt_edge_coords[single_crossing_idx, :] - refined_edge_coords[single_crossing_idx, :])**2, axis=1))
mdist = np.mean(dist)
return opt_xfm
return opt_xfm, mdist
def image_props(img, img_resolution):
......@@ -455,7 +514,7 @@ def justify_bbox(bbox, aspect_ratio, justify):
def initialise_xfm(image, image_resolution, contour, tol=0.05, justify=None):
'''Calculate image and contour bounding boxes, then estimate the transform (xfm) based on bounding box corners.'''
"""Calculate image and contour bounding boxes, then estimate the transform (xfm) based on bounding box corners."""
brainmask = segment_foreground(image)
......@@ -468,20 +527,20 @@ def initialise_xfm(image, image_resolution, contour, tol=0.05, justify=None):
bbox = justify_bbox(image_p["bbox"], contour_p["aspect_ratio"], justify)
# exclude voxels from brainmask that are outside the bounding box
# # exclude voxels from brainmask that are outside the bounding box
bbox = np.round(np.array(bbox) / image_resolution).astype(int)
# bbox = np.round(np.array(bbox) / image_resolution).astype(int)
cropped_brainmask = brainmask.copy()
cropped_brainmask[:, : bbox[1]] = 0
cropped_brainmask[:, bbox[3] :] = 0
# cropped_brainmask = brainmask.copy()
# cropped_brainmask[:, : bbox[1]] = 0
# cropped_brainmask[:, bbox[3] :] = 0
image_p = image_props(cropped_brainmask, image_resolution)
# image_p = image_props(cropped_brainmask, image_resolution)
# recalculate justified bounding box on cropped brainmask
# force image bbox aspect ratio to be the same as contour bbox, and justify left or right.
# # recalculate justified bounding box on cropped brainmask
# # force image bbox aspect ratio to be the same as contour bbox, and justify left or right.
bbox = justify_bbox(image_p["bbox"], contour_p["aspect_ratio"], justify)
# bbox = justify_bbox(image_p["bbox"], contour_p["aspect_ratio"], justify)
centroid = (
(bbox[0] + bbox[2]) / 2,
......@@ -497,7 +556,9 @@ def initialise_xfm(image, image_resolution, contour, tol=0.05, justify=None):
}
if np.abs(image_p["aspect_ratio"] - contour_p["aspect_ratio"]) > tol:
raise RuntimeError("Slide and chart aspect ratios appear to be different. Consider correcting aspect-ratio with left/right justification.")
raise RuntimeError(
"Slide and chart aspect ratios appear to be different. Consider correcting aspect-ratio with left/right justification."
)
idx = np.array([[1, 2], [1, 0], [3, 2], [3, 0]])
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment