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

Added support for extended chart to chart_reg

parent 3ced85c3
......@@ -20,7 +20,6 @@ import json
import yaml
from scipy.ndimage.interpolation import affine_transform
from slider.external import neurolucida
from slider import util
import glymur
......@@ -62,50 +61,25 @@ def estimate_field(img):
return 'dark'
def enhance(img0, kernel_size=None, lower_percentile=2, upper_percentile=98, sigma=5):
'''Enhance image exposure and contrast'''
if kernel_size is None:
h, w = img0.shape
kernel_size=(h//8, w//8)
img0 = exposure.equalize_adapthist(img0, kernel_size=kernel_size)
p_lower, p_upper = np.percentile(img0, (lower_percentile, upper_percentile))
img0 = exposure.rescale_intensity(img0, in_range=(p_lower, p_upper))
class Slide:
img0 = filters.gaussian(img0, sigma)
return img0
def segment_foreground(img, marker_threshold=(0.02, 0.2), min_component_size=10000, selem=morphology.disk(10)):
'''Segment image foreground from background'''
contours = None
cells = None
# calculate elevation map
elevation_map = filters.sobel(img)
# extract background and foreground features
markers = np.zeros_like(img)
markers[img < marker_threshold[0]] = 1
markers[img > marker_threshold[1]] = 2
def __init__(self, image, resolution) -> None:
self.data = image
self.resolution = resolution
# segment using watershed algorithm
labels = segmentation.watershed(elevation_map, markers)
labels = ndi.binary_fill_holes(labels-1)+1
@property
def shape(self):
return self.data.shape
# find connected components and exclude components < min_component_size
cc = measure.label(labels, background=1)
p = measure.regionprops(cc, img)
@classmethod
def from_jp2(cls, fname, resolution):
im = glymur.Jp2k(fname)
return cls(im, resolution)
ridx = np.where([p0.area>min_component_size for p0 in p])[0]
brainmask = np.sum(np.stack([(cc == p[r].label) for r in ridx], axis=2), axis=2)
# erode and dilate to clean up edges
brainmask = morphology.binary_erosion(brainmask, selem=selem)
brainmask = morphology.binary_dilation(brainmask, selem=selem)
return brainmask
def register_chart_to_slide(chart, slide, slide_res, outdir, boundary_key=None, config=None, do_plots=None):
......@@ -130,16 +104,13 @@ def register_chart_to_slide(chart, slide, slide_res, outdir, boundary_key=None,
# load chart
contour, cells = neurolucida.read(chart)
boundary_idx = np.where([f.lower()==boundary_key.lower() for f,_ in contour])[0]
# if boundary_idx.size > 1:
# raise RuntimeError(f'Repeated entries of boundary key in chart: {boundary_key}')
chart = util.Chart.from_neurolucida_ascii(chart)
# edge_crds = contour[boundary_idx[0]][1][:, :2] * [1, -1]
outline = chart.get_contours(name=boundary_key, closed=True)
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}')
edge_crds = [contour[bidx][1][:, :2] * [1, -1] for bidx in boundary_idx]
edge_crds = np.concatenate([e for e in edge_crds[:2]], axis=0)
edge_crds = np.concatenate([x.points[:, :2] for x in outline]) * [1, -1]
# load slide, convert to grayscale, and invert if light-field
......@@ -189,7 +160,7 @@ def register_chart_to_slide(chart, slide, slide_res, outdir, boundary_key=None,
# apply opt-xfm to contours and cells and save
contour_xfm = [
(k, apply_xfm(init, v[:, :2] * [1, -1]).tolist()) for k, v in contour
(contour.name, apply_xfm(init, contour.points[:, :2] * [1, -1]).tolist(), contour.closed) for contour in chart.contours
]
if do_plots:
......@@ -198,7 +169,9 @@ def register_chart_to_slide(chart, slide, slide_res, outdir, boundary_key=None,
extent = np.array([0, img.shape[1], img.shape[0], 0]) * slide_res
ax.imshow(img, extent=extent, cmap='gray')
for name, coords in contour_xfm:
for name, coords, closed in contour_xfm:
if not closed: continue
coords = np.array(coords)
......@@ -227,10 +200,56 @@ def register_chart_to_slide(chart, slide, slide_res, outdir, boundary_key=None,
with open(f"{outdir}/contour.json", "w") as fp:
json.dump(contour_xfm, fp, indent=4)
if cells.size > 0:
cells_xfm = apply_xfm(init, cells[:, :2] * [1, -1])
with open(f"{outdir}/cells.json", "w") as fp:
json.dump(cells_xfm.tolist(), fp, indent=4)
# if cells.size > 0:
# cells_xfm = apply_xfm(init, cells[:, :2] * [1, -1])
# with open(f"{outdir}/cells.json", "w") as fp:
# json.dump(cells_xfm.tolist(), fp, indent=4)
def enhance(img0, kernel_size=None, lower_percentile=2, upper_percentile=98, sigma=5):
'''Enhance image exposure and contrast'''
if kernel_size is None:
h, w = img0.shape
kernel_size=(h//8, w//8)
img0 = exposure.equalize_adapthist(img0, kernel_size=kernel_size)
p_lower, p_upper = np.percentile(img0, (lower_percentile, upper_percentile))
img0 = exposure.rescale_intensity(img0, in_range=(p_lower, p_upper))
img0 = filters.gaussian(img0, sigma)
return img0
def segment_foreground(img, marker_threshold=(0.02, 0.2), min_component_size=10000, selem=morphology.disk(30)):
'''Segment image foreground from background'''
# calculate elevation map
elevation_map = filters.sobel(img)
# extract background and foreground features
markers = np.zeros_like(img)
markers[img < marker_threshold[0]] = 1
markers[img > marker_threshold[1]] = 2
# segment using watershed algorithm
labels = segmentation.watershed(elevation_map, markers)
labels = ndi.binary_fill_holes(labels-1)+1
# find connected components and exclude components < min_component_size
cc = measure.label(labels, background=1)
p = measure.regionprops(cc, img)
ridx = np.where([p0.area>min_component_size for p0 in p])[0]
brainmask = np.sum(np.stack([(cc == p[r].label) for r in ridx], axis=2), axis=2)
# erode and dilate to clean up edges
brainmask = morphology.binary_erosion(brainmask, selem=selem)
brainmask = morphology.binary_dilation(brainmask, selem=selem)
return brainmask
def refine_edge_coord(img, img_res, edge_coords, normals, do_plots=True):
......@@ -411,7 +430,7 @@ def point_props(pnts, do_plots=False, plot_name=None):
def init_scale(img, img_resolution, crd):
def init_scale(img, img_resolution, crd, tol=0.05):
img_p = image_props(img, img_resolution, do_plots=True)
crd_p = point_props(crd, do_plots=True)
......@@ -446,7 +465,7 @@ def init_scale(img, img_resolution, crd):
# print(img_p)
# print(crd_p)
if np.abs(img_p['aspect_ratio'] - crd_p['aspect_ratio']) > 0.005:
if np.abs(img_p['aspect_ratio'] - crd_p['aspect_ratio']) > tol:
raise RuntimeError('Slide and chart aspect ratios appear to be different.')
idx = np.array([[1, 2], [1, 0], [3, 2], [3, 0]])
......
......@@ -16,6 +16,10 @@
#
import os.path as op
from dataclasses import dataclass
from typing import Optional, List
import numpy as np
from slider.external import neurolucida
def get_slider_dir() -> str:
rpath = op.realpath(op.dirname(op.abspath(__file__)))
......@@ -30,3 +34,61 @@ def get_resource(name) -> str:
r = op.join(get_resource_path(), name)
# asrt.assert_file_exists(r)
return r
@dataclass
class ChartContour:
'''Container for chart contour'''
name: str
closed: bool
color: str
resolution: float
points: np.array
guid: Optional[str] = None
def __repr__(self):
r, c = self.points.shape
return f"ChartContour(name='{self.name}', guid={self.guid}, closed={self.closed}, color='{self.color}', resolution={self.resolution}, points=[{r}x{c} np.array])"
@dataclass
class ChartCell:
'''Container for chart cell'''
color: str
point: np.array
@dataclass
class Chart:
'''Container for chart'''
contours: List[ChartContour] = None
cells: List[ChartCell] = None
@property
def n_contours(self):
return len(self.contours) if self.contours is not None else 0
@property
def n_cells(self):
return len(self.cells) if self.cells is not None else 0
def get_contours(self, name=None, closed=None):
cnt = self.contours
if name is not None: cnt = [c for c in cnt if c.name==name]
if closed is not None: cnt = [c for c in cnt if c.closed==closed]
return cnt
def get_contour_names(self):
return [contour.name for contour in c.contours] if self.n_contours > 0 else None
@classmethod
def from_neurolucida_ascii(cls, fname):
contours, cells = neurolucida.read(fname)
return cls(
[ChartContour(**cnt) for cnt in contours],
[ChartCell(**cell) for cell in cells]
)
def __repr__(self):
return f"Chart(contours=[{self.n_contours}x ChartContour], cells=[{self.n_cells}x ChartCell])"
\ No newline at end of file
Markdown is supported
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