chart_reg.py 6.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#!/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

17
def register_chart_to_slide(chart, slide, slide_res, out, rlevel=4, boundary_key='outline'):
18
19
20
21

    # load chart
    contour, cells = neurolucida.read(chart)

22
23
    slide_res = slide_res * (2**rlevel)

24
25
26
27
28
29
    # 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]
30
    init = _init_scale(img, slide_res, edge_crds, verbose=True, plots=False)
31
32
33
34
35
36
37
38
    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)

Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
39
    # refine edge_coords (to image)
40
    refined_edge_coords = _refine_edge_coord(img, slide_res, xfm_edge_coords, init_nrmls)
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

    # 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

57
58
59
    contour_xfm = {
        k: _apply_xfm(opt, v[:, :2] * [1, -1]).tolist() for k, v in contour.items()
    }
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
60
61
62
63
64
65
66
67
68
    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)


69
def _refine_edge_coord(img, img_res, edge_coords, normals):
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
70
71
72
73
74
75
76
    """
    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
77

78
    line_smpl = np.arange(-20, 20) * img_res
79

80
81
    line_x = edge_x[:, np.newaxis] + nrml_x[:, np.newaxis] * line_smpl
    line_y = edge_y[:, np.newaxis] + nrml_y[:, np.newaxis] * line_smpl
82

83
84
85
    line_x = np.round(line_x / img_res).astype(int)
    line_y = np.round(line_y / img_res).astype(int)
    
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

    # 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)

117
    return refined_edge_coords * img_res
118
119
120



121
def _img_props(img, img_resolution, verbose=False, plots=False):
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144

    # 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]

145
    bbox = np.array(p.bbox) * img_resolution
146
147
148
149
150
151
152

    centroid = (
        (bbox[0] + bbox[2])/2,
        (bbox[1] + bbox[3])/2,
    )

    # TODO: make plots work again
153
    # if plots:
154

155
156
    #     plt.imshow(cc == p.label)
    #     plt.colorbar()
157

158
159
160
161
162
163
164
165
166
    #     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)
167

168
169
    #     plt.plot(centroid[1], centroid[0], 'ro')
    #     plt.show()
170

171
172
173
174
175
    return {
        'bbox': bbox, 
        'bbox_centroid': centroid, 
        'shape': (bbox[2]-bbox[0], bbox[3]-bbox[1])
        }
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219


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])}


220
def _init_scale(img, img_resolution, crd, verbose=False, plots=False):
221

222
    img_p = _img_props(img, img_resolution)
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
    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