Skip to content
GitLab
Commit 607a6b4f authored by Sean Fitzgibbon's avatar Sean Fitzgibbon
Browse files

Improved plotting

parent e2a6ee31
Hide whitespace changes
Inline Side-by-side
slider/chart_reg.py
View file @ 607a6b4f
......@@ -34,7 +34,6 @@ from skimage import (
transform,
)
from typing import Optional
from slider import util
......@@ -43,14 +42,14 @@ def register_chart_to_slide(
slide: str,
slide_res: float,
outdir: str,
boundary_key: Optional[str] = None,
boundary_key: Optional[str] = ("Outline", "outline", "Section Outline"),
config: Optional[str] = None,
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.
......@@ -60,7 +59,7 @@ def register_chart_to_slide(
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")
......@@ -84,11 +83,18 @@ def register_chart_to_slide(
chart = util.Chart.from_neurolucida_ascii(chart)
outline = chart.get_contours(name=boundary_key, closed=True)
if isinstance(boundary_key, str):
outline = chart.get_contours(name=boundary_key, closed=True)
else:
outline = [chart.get_contours(name=k, closed=True) for k in boundary_key]
outline = [item for sublist in outline for item in sublist]
if len(outline) < 1:
raise ValueError(f"Boundary key {boundary_key} not found in chart")
if len(outline) > 1:
warnings.warn(f"Repeated entries of boundary key in chart: {boundary_key}")
warnings.warn(
f"Repeated entries of boundary key in chart ({boundary_key}). They will be merged."
)
edge_crds = [x.points[:, :2] * [1, -1] for x in outline]
edge_crds_cat = np.concatenate(edge_crds)
......@@ -123,19 +129,18 @@ def register_chart_to_slide(
np.savetxt(f"{outdir}/chart-to-image-init.xfm", init_xfm.params)
# refine alignment (to mask edges)
mask = img_props['mask']
opt_xfm, mdist = optimise_xfm(mask, slide_res, edge_crds_cat[::3,:], init_xfm)
# optimise alignment
# The `optimise_xfm` function takes a mask and a set of coordinates, and returns the optimal
# transformation to map the coordinates to the mask.
mask = img_props["mask"]
opt_xfm, mdist, iter_props = 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
......@@ -168,10 +173,7 @@ 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))
fig, ax = plt.subplots(4, 2, figsize=(20, 30))
ax = np.ravel(ax)
# chart bounding box
......@@ -194,33 +196,101 @@ def register_chart_to_slide(
)
plot_box(ax[1], img_props["bbox"])
# normals (first iteration)
plot_slide(
ax[2], img_props["mask"], slide_res, "Brainmask + normals: first iteration"
)
line_x, line_y = iter_props[0]["normals"]
contour_x, contour_y = iter_props[0]["contour_init"]
single_crossing_idx = iter_props[0]["single_crossing_idx"]
plot_normals(
ax[2],
line_x[~single_crossing_idx, :],
line_y[~single_crossing_idx, :],
contour_x,
contour_y,
color="y",
)
plot_normals(
ax[2],
line_x[single_crossing_idx, :],
line_y[single_crossing_idx, :],
contour_x,
contour_y,
color="b",
)
# normals (last iteration)
plot_slide(
ax[3], img_props["mask"], slide_res, "Brainmask + normals: last iteration"
)
line_x, line_y = iter_props[-1]["normals"]
contour_x, contour_y = iter_props[-1]["contour_init"]
single_crossing_idx = iter_props[-1]["single_crossing_idx"]
plot_normals(
ax[3],
line_x[~single_crossing_idx, :],
line_y[~single_crossing_idx, :],
contour_x,
contour_y,
color="y",
)
plot_normals(
ax[3],
line_x[single_crossing_idx, :],
line_y[single_crossing_idx, :],
contour_x,
contour_y,
color="b",
)
# distance plot
ax[4].plot(mdist)
ax[4].set_xlabel("iteration")
ax[4].set_ylabel("mean residual distance")
ax[5].set_axis_off()
# boundary plots
plot_slide(ax[2], img, slide_res, title="Boundary alignment")
plot_slide(ax[6], img, slide_res, title="Boundary alignment")
for contour in edge_crds:
plot_contour(
ax[2], apply_xfm(init_xfm, contour), color=(1, 0, 0), marker="."
ax[6], apply_xfm(init_xfm, contour), color=(1, 0, 0), marker="."
)
plot_contour(
ax[2], apply_xfm(opt_xfm, contour), color=(0, 1, 0), marker="."
ax[6], apply_xfm(opt_xfm, contour), color=(0, 1, 0), marker="."
)
ax[2].legend(["boundary_init", "boundary_optimised"])
ax[6].legend(["boundary_init", "boundary_optimised"])
# aligned chart
plot_slide(ax[3], img, slide_res, title="Aligned Chart")
plot_slide(ax[7], img, slide_res, title="Aligned Chart")
cmap = plt.get_cmap("tab10")
for idx, (name, coords, closed) in enumerate(contour_xfm):
if not closed:
continue
plot_contour(ax[3], coords, name, color=cmap(idx))
plot_contour(ax[7], coords, name, color=cmap(idx))
fig.savefig(f"{outdir}/alignment.png", bbox_inches="tight", dpi=300)
def plot_normals(ax, line_x, line_y, contour_x, contour_y, color="b"):
for line_x0, line_y0 in zip(line_x, line_y):
ax.plot(line_x0, line_y0, ".", color=color)
ax.plot(contour_x, contour_y, "r.")
def plot_box(ax, bbox, edgecolor="r", facecolor="none", linewidth=1):
rect = patches.Rectangle(
(bbox[1], bbox[0]),
......@@ -338,16 +408,20 @@ def segment_foreground(
def optimise_xfm(mask, mask_resolution, contour, xfm, n_iter=100):
mdist = np.zeros(n_iter)
iter_props = [None] * n_iter
# TODO: move line extent values to config
max_line_extent = 1.5
min_line_extent = 0.01
min_line_extent = 0.05
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
xfm, mdist[idx], iter_props[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, iter_props
def optimise_xfm_worker(image, image_resolution, contour, xfm_init, line_extent=1.5):
......@@ -383,44 +457,7 @@ def optimise_xfm_worker(image, image_resolution, contour, xfm_init, line_extent=
line_int = brainmask[np.clip(line_y, 0, y - 1), np.clip(line_x, 0, x - 1)]
# 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))
extent = np.array([0, image.shape[1], image.shape[0], 0])
ax.imshow(brainmask, extent=extent, cmap="gray")
for line_x0, line_y0 in zip(line_x[single_crossing_idx, :], line_y[single_crossing_idx, :]):
ax.plot(line_x0, line_y0, "b.")
ax.plot(edge_init_x, edge_init_y, "r.")
ax.set_xlabel("mm")
ax.set_ylabel("mm")
ax.set_title("Brainmask + normals")
plt.show()
# exlude constant rows
# constant_idx = np.all(line_int == line_int[:, 0][:, np.newaxis], axis=1)
single_crossing_idx = np.sum(np.abs(np.diff(line_int, axis=1)), axis=1) == 1
min_idx = np.argmax(np.abs(np.diff(line_int, axis=-1)), axis=-1)
......@@ -433,16 +470,33 @@ def optimise_xfm_worker(image, image_resolution, contour, xfm_init, line_extent=
)
opt_xfm = transform.SimilarityTransform()
opt_xfm.estimate(contour[single_crossing_idx, :], refined_edge_coords[single_crossing_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))
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, mdist
iter_props = {
"single_crossing_idx": single_crossing_idx,
"normals": (line_x * image_resolution, line_y * image_resolution),
"contour_init": (edge_init_x, edge_init_y),
}
return opt_xfm, mdist, iter_props
def image_props(img, img_resolution):
......@@ -501,16 +555,18 @@ def justify_bbox(bbox, aspect_ratio, justify):
new_w = w * scale_factor
new_bbox = bbox.copy()
if justify == "right":
bbox[1] += w - new_w
new_bbox[1] += w - new_w
elif justify == "left":
bbox[3] -= w - new_w
new_bbox[3] -= w - new_w
else:
raise ValueError(
f'Unknown value for justify ("{justify}"). Can only be "left" or "right".'
)
return bbox
return new_bbox
def initialise_xfm(image, image_resolution, contour, tol=0.05, justify=None):
......@@ -521,6 +577,49 @@ def initialise_xfm(image, image_resolution, contour, tol=0.05, justify=None):
image_p = image_props(brainmask, image_resolution)
contour_p = point_props(contour)
# If the aspect-ratio of the image and contour bounding boxes are not equal it typically means that the image is bilateral whilst the contour is lateralised.
# Thust the contour bbox should left or right justified within the image bbox.
def aspect_is_equal(x, y):
return np.abs(x["aspect_ratio"] - y["aspect_ratio"]) < tol
if (not aspect_is_equal(image_p, contour_p)) & (justify is None):
# The justification is not set by the user, will attempt left and right and see which has best fit.
warnings.warn(
"Slide and chart aspect ratios appear to be different, and direction of justification is not set. Will attempt to estimate direction."
)
bbox_idx = np.array([[1, 2], [1, 0], [3, 2], [3, 0]])
src = np.array(contour_p["bbox"])[bbox_idx]
goodness_fit = np.zeros(2)
directions = ("left", "right")
for dir_idx, dir in enumerate(directions):
bbox = justify_bbox(image_p["bbox"], contour_p["aspect_ratio"], dir)
dest = np.array(bbox)[bbox_idx]
xfm = transform.SimilarityTransform()
xfm.estimate(src, dest)
# calc goodness of fit (number of single-edge-crossing normals)
# fit is determined by the number of contour normals that cross the image mask boundary only once
single_crossing_idx = optimise_xfm_worker(
brainmask, image_resolution, contour, xfm, line_extent=1
)[2]["single_crossing_idx"]
goodness_fit[dir_idx] = np.sum(single_crossing_idx)
# maximise goodness-of-fit
justify = directions[np.argmax(goodness_fit)]
print(
f"Estimated that the contour should be {justify} justified to the image: {goodness_fit}."
)
if justify is not None:
# force image bbox aspect ratio to be the same as contour bbox, and justify left or right.
......@@ -555,7 +654,7 @@ def initialise_xfm(image, image_resolution, contour, tol=0.05, justify=None):
"mask": brainmask,
}
if np.abs(image_p["aspect_ratio"] - contour_p["aspect_ratio"]) > tol:
if not aspect_is_equal(image_p, contour_p):
raise RuntimeError(
"Slide and chart aspect ratios appear to be different. Consider correcting aspect-ratio with left/right justification."
)
......
slider/util.py
View file @ 607a6b4f
......@@ -21,6 +21,11 @@ from typing import Optional, List
import numpy as np
from slider.external import neurolucida
import glymur
import collections
import six
def is_nonstring_iterable(arg):
return (isinstance(arg, collections.Iterable) and not isinstance(arg, six.string_types))
def get_slider_dir() -> str:
rpath = op.realpath(op.dirname(op.abspath(__file__)))
......
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