#!/usr/bin/env python import os import os.path as op import json from scipy.ndimage.interpolation import affine_transform from slider.external import neurolucida import glymur from skimage.measure import regionprops, label from skimage.color import rgb2gray from skimage import filters, transform from matplotlib import patches import numpy as np import matplotlib.pyplot as plt def register_chart_to_slide(chart, slide, slide_res, out, rlevel=4, boundary_key='outline'): # load chart contour, cells = neurolucida.read(chart) slide_res = slide_res * (2**rlevel) # load slide jp2 = glymur.Jp2k(slide) img = jp2.read(rlevel=rlevel) # initial scaling based on boundng boxes edge_crds = contour[boundary_key][:, :2] * [1, -1] init = _init_scale(img, slide_res, edge_crds, verbose=True, plots=False) print(init) print( f'Rotation: {init.rotation}, Translation: {init.translation}, Scale: {init.scale}, Shear: {init.shear}') # calculate normal line to boundary points xfm_edge_coords = _apply_xfm(init, edge_crds) init_nrmls = normal(xfm_edge_coords, plot=False, distance=1, verbose=True) # refine edge_coords (to image) refined_edge_coords = _refine_edge_coord(img, slide_res, xfm_edge_coords, init_nrmls) # estimate opimised affine transform opt = transform.AffineTransform() opt.estimate(edge_crds, refined_edge_coords) print(opt) print( f'Rotation:\t{opt.rotation}\nTranslation:\t{opt.translation}\nScale:\t\t{list(opt.scale)}\nShear:\t\t{opt.shear}') # save opt transform if not op.exists(out): os.makedirs(out) np.savetxt(f'{out}/chart-to-image.xfm', opt.params) # apply opt-xfm to contours and cells and save contour_xfm = { k: _apply_xfm(opt, v[:, :2] * [1, -1]).tolist() for k, v in contour.items() } with open(f'{out}/contour.json', 'w') as fp: json.dump(contour_xfm, fp, indent=4) if cells.size > 0: cells_xfm = _apply_xfm(cells[:, :2] * [1, -1]) with open(f'{out}/cells.json', 'w') as fp: json.dump(cells_xfm, fp, indent=4) def _refine_edge_coord(img, img_res, edge_coords, normals): """ Refine edge_coord by sampling binarised image along normal (to edge) and looking for big step change. """ # calculate normal line (to edge_coords) edge_x, edge_y = np.round(edge_coords.T).astype(int) nrml_x, nrml_y = normals.T line_smpl = np.arange(-20, 20) * img_res line_x = edge_x[:, np.newaxis] + nrml_x[:, np.newaxis] * line_smpl line_y = edge_y[:, np.newaxis] + nrml_y[:, np.newaxis] * line_smpl line_x = np.round(line_x / img_res).astype(int) line_y = np.round(line_y / img_res).astype(int) # find threshold for background image = rgb2gray(img) threshold_value = filters.threshold_otsu(image) print(f'threshold={threshold_value}') # binarise foreground # NOTE: this assumes a bright background! labeled_foreground = (image < threshold_value).astype(int) # calculate connected components cc = label(labeled_foreground, background=0) # get region properties for each cc and select the largest cc (in terms of area) p = regionprops(cc, image) ridx = np.argmax([p0.area for p0 in p]) p = p[ridx] cc0 = (cc == p.label) # plt.imshow(cc0) line_int = cc0[line_y, line_x] min_idx = np.argmax(np.diff(line_int, axis=-1), axis=-1) refined_edge_coords = np.stack([ line_x[np.arange(min_idx.size), min_idx], line_y[np.arange(min_idx.size), min_idx], ], axis=-1) return refined_edge_coords * img_res def _img_props(img, img_resolution, verbose=False, plots=False): # convert to grayscale image = rgb2gray(img) # find threshold for background threshold_value = filters.threshold_otsu(image) if verbose: print(f'threshold={threshold_value}') # binarise foreground # NOTE: this assumes a bright background! labeled_foreground = (image < threshold_value).astype(int) # calculate connected components cc = label(labeled_foreground, background=0) # get region properties for each cc and select the largest cc (in terms of area) p = regionprops(cc, image) ridx = np.argmax([p0.area for p0 in p]) p = p[ridx] bbox = np.array(p.bbox) * img_resolution centroid = ( (bbox[0] + bbox[2])/2, (bbox[1] + bbox[3])/2, ) # TODO: make plots work again # if plots: # plt.imshow(cc == p.label) # plt.colorbar() # rect = patches.Rectangle( # (bbox[1], bbox[0]), # bbox[3]-bbox[1], # bbox[2]-bbox[0], # linewidth=1, # edgecolor='r', # facecolor='none' # ) # plt.gca().add_patch(rect) # plt.plot(centroid[1], centroid[0], 'ro') # plt.show() return { 'bbox': bbox, 'bbox_centroid': centroid, 'shape': (bbox[2]-bbox[0], bbox[3]-bbox[1]) } def _pnt_props(pnts, verbose=False, plots=False): x = pnts[:, 0] y = pnts[:, 1] bbox = ( np.min(y), np.min(x), np.max(y), np.max(x) ) centroid = ( (bbox[0] + bbox[2])/2, (bbox[1] + bbox[3])/2, ) if verbose: print(f'bbox={bbox}') print(f'centroid={centroid}') if plots: plt.plot(x, y, 'b.') plt.axis('equal') plt.gca().invert_yaxis() rect = patches.Rectangle( (bbox[1], bbox[0]), bbox[3]-bbox[1], bbox[2]-bbox[0], linewidth=1, edgecolor='r', facecolor='none' ) plt.gca().add_patch(rect) plt.plot(centroid[1], centroid[0], 'ro') return {'bbox': bbox, 'bbox_centroid': centroid, 'shape': (bbox[2]-bbox[0], bbox[3]-bbox[1])} def _init_scale(img, img_resolution, crd, verbose=False, plots=False): img_p = _img_props(img, img_resolution) crd_p = _pnt_props(crd) if verbose: print(f'image props = {img_p}') print(f'chart props = {crd_p}') idx = np.array([[1, 2], [1, 0], [3, 2], [3, 0]]) src = np.array(crd_p['bbox'])[idx] dest = np.array(img_p['bbox'])[idx] a = transform.AffineTransform() a.estimate(src, dest) return a def _apply_xfm(xfm, pnts): return (xfm.params @ np.concatenate((pnts, np.ones((pnts.shape[0], 1))), axis=-1).T).T[:, :2] def normal(pnts, distance=100, plot=False, verbose=False): d_pnts = np.roll(pnts, 1, axis=0) - np.roll(pnts, -1, axis=0) normal = np.stack([d_pnts[:, 1], -d_pnts[:, 0]], axis=-1) normal = normal / np.linalg.norm(normal, axis=-1, keepdims=True) return normal