chart_reg.py 10.7 KB
Newer Older
1
#!/usr/bin/env python
2
#
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
3
# Copyright 2021 Sean Fitzgibbon, University of Oxford
4
#
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
5
6
7
#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at
8
#
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
9
#        http://www.apache.org/licenses/LICENSE-2.0
10
#
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
11
12
13
14
15
#    Unless required by applicable law or agreed to in writing, software
#    distributed under the License is distributed on an "AS IS" BASIS,
#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#    See the License for the specific language governing permissions and
#    limitations under the License.
16
#
17
18
19
import os
import os.path as op
import json
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
20
import yaml
21
22
23

from scipy.ndimage.interpolation import affine_transform
from slider.external import neurolucida
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
24
from slider import util
25
26
27
import glymur

from skimage.measure import regionprops, label
28
29
from skimage.color import rgb2gray, rgb2hsv
from skimage import filters, transform, exposure
30
31
32
33
from matplotlib import patches
import numpy as np
import matplotlib.pyplot as plt

34
35
36
DO_PLOTS = False
OUTDIR = None

37

Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
38
def register_chart_to_slide(chart, slide, slide_res, out, boundary_key=None, config=None):
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
39
40

    if config is None:
41
        config = util.get_resource("chart.yaml")
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
42

43
    with open(config, "r") as f:
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
44
        config = yaml.load(f, Loader=yaml.FullLoader)
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
45

46
    rlevel = config["slide"]["resolution_level"]
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
47
48
49
50
51

    if boundary_key is None:
        boundary_key = config["chart"]["boundary_key"]

    print(boundary_key)
52

53
    # set global variables required for summary plots
54
    plots = config["general"]["plots"]
55
56
57

    DO_PLOTS = plots
    OUTDIR = out
58
59
60
61

    #  create output dir
    if not op.exists(out):
        os.makedirs(out)
62

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

66
    slide_res = slide_res * (2 ** rlevel)
67

68
69
70
71
72
    # load slide
    jp2 = glymur.Jp2k(slide)
    img = jp2.read(rlevel=rlevel)

    # initial scaling based on boundng boxes
73
74
75
76
77
78
    boundary_idx = np.where([f==boundary_key for f,_ in contour])[0]
    if boundary_idx.size > 1:
        raise RuntimeError(f'Repeated entries of boundary key in chart: {boundary_key}')
    print(boundary_idx)

    edge_crds = contour[boundary_idx[0]][1][:, :2] * [1, -1]
79
    init = _init_scale(img, slide_res, edge_crds, verbose=True)
80
    print(init)
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
81
82
83
    # print(
    #     f"Rotation: {init.rotation}, Translation: {init.translation}, Scale: {init.scale}, Shear: {init.shear}"
    # )
84
85
86

    # calculate normal line to boundary points
    xfm_edge_coords = _apply_xfm(init, edge_crds)
87
    init_nrmls = normal(xfm_edge_coords)
88

Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
89
    # refine edge_coords (to image)
90
91
92
    refined_edge_coords = _refine_edge_coord(
        img, slide_res, xfm_edge_coords, init_nrmls
    )
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
93
94

    # estimate opimised affine transform
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
95
    opt = transform.SimilarityTransform()
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
96
97
    opt.estimate(edge_crds, refined_edge_coords)
    print(opt)
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
98
99
100
    # print(
    #     f"Rotation:\t{opt.rotation}\nTranslation:\t{opt.translation}\nScale:\t\t{list(opt.scale)}\nShear:\t\t{opt.shear}"
    # )
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
101
102

    # save opt transform
103
    np.savetxt(f"{out}/chart-to-image.xfm", opt.params)
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
104
105
106

    # apply opt-xfm to contours and cells and save

107
108
109
    contour_xfm = [
        (k, _apply_xfm(opt, v[:, :2] * [1, -1]).tolist()) for k, v in contour
    ]
110
111
112
113
114
115
116

    if DO_PLOTS:
        fig, ax = plt.subplots(figsize=(10, 10))

        extent = np.array([0, img.shape[1], img.shape[0], 0]) * slide_res
        ax.imshow(img, extent=extent)

117
        for name, coords in contour_xfm:
118
119
120

            coords = np.array(coords)

121
122
123
124
125
126
127
            cog = np.mean(coords, axis=0)
            ax.text(cog[0], cog[1], name)
            ax.plot(coords[:, 0], coords[:, 1], "k-")

        ax.set_xlabel("mm")
        ax.set_ylabel("mm")
        ax.set_title("Aligned Chart")
128

129
        fig.savefig(f"{OUTDIR}/aligned_chart.png", bbox_inches="tight", dpi=300)
130

131
        # fig, ax = plt.subplots(1, 5, figsize=(25, 5))
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
132
133
134
135
136
137
138
139
140
141
142
        # ax[0].imshow(plt.imread(f'{OUTDIR}/chart_bounding_box.png'))
        # ax[1].imshow(plt.imread(f'{OUTDIR}/image_bounding_box.png'))
        # ax[2].imshow(plt.imread(f'{OUTDIR}/normals.png'))
        # ax[3].imshow(plt.imread(f'{OUTDIR}/refined_coords.png'))
        # ax[4].imshow(plt.imread(f'{OUTDIR}/aligned_chart.png'))

        # for a in ax:
        #     a.axis('off')

        # fig.savefig(f'{OUTDIR}/alignment.png', bbox_inches='tight', dpi=150)

143
    with open(f"{out}/contour.json", "w") as fp:
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
144
145
146
        json.dump(contour_xfm, fp, indent=4)

    if cells.size > 0:
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
147
        cells_xfm = _apply_xfm(opt, cells[:, :2] * [1, -1])
148
        with open(f"{out}/cells.json", "w") as fp:
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
149
            json.dump(cells_xfm.tolist(), fp, indent=4)
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
150
151


152
def _refine_edge_coord(img, img_res, edge_coords, normals):
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
153
154
155
156
157
    """
    Refine edge_coord by sampling binarised image along normal (to edge) and looking for big step change.
    """

    # calculate normal line (to edge_coords)
158
    edge_x, edge_y = edge_coords.T
Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
159
    nrml_x, nrml_y = normals.T
160

161
    # TODO: move these line extents to the config file
162
    line_smpl = np.linspace(-0.03, 0.15, 20)
163

164
165
    line_x = edge_x[:, np.newaxis] + nrml_x[:, np.newaxis] * line_smpl
    line_y = edge_y[:, np.newaxis] + nrml_y[:, np.newaxis] * line_smpl
166

167
168
169
170
    # find threshold for background
    image = rgb2gray(img)
    threshold_value = filters.threshold_otsu(image)

171
    print(f"threshold={threshold_value}")
172
173
174
175
176
177
178
179
180
181
182
183
184

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

185
    cc0 = cc == p.label
186
187
188
189
190
191
192
193

    if DO_PLOTS:
        fig, ax = plt.subplots(figsize=(20, 30))

        extent = np.array([0, img.shape[1], img.shape[0], 0]) * img_res
        ax.imshow(cc0, extent=extent)

        for line_x0, line_y0 in zip(line_x, line_y):
194
            ax.plot(line_x0, line_y0, "b.")
195

196
        ax.plot(edge_x, edge_y, "r.")
197

198
199
200
        ax.set_xlabel("mm")
        ax.set_ylabel("mm")
        ax.set_title("Largest connected component + normals")
201

202
        fig.savefig(f"{OUTDIR}/normals.png", bbox_inches="tight", dpi=300)
203
204
205
206
207

    # sample image along normal line

    line_x = np.round(line_x / img_res).astype(int)
    line_y = np.round(line_y / img_res).astype(int)
208
209
210
211
212

    line_int = cc0[line_y, line_x]

    min_idx = np.argmax(np.diff(line_int, axis=-1), axis=-1)

213
214
215
216
217
218
219
    refined_edge_coords = np.stack(
        [
            line_x[np.arange(min_idx.size), min_idx] * img_res,
            line_y[np.arange(min_idx.size), min_idx] * img_res,
        ],
        axis=-1,
    )
220

221
222
223
224
225
226
227
    # TODO: plot edge_coords + refined_edge_coords
    if DO_PLOTS:
        fig, ax = plt.subplots(figsize=(20, 30))

        extent = np.array([0, img.shape[1], img.shape[0], 0]) * img_res
        ax.imshow(img, extent=extent)

228
229
        ax.plot(edge_x, edge_y, "r.-")
        ax.plot(refined_edge_coords[:, 0], refined_edge_coords[:, 1], "g.-")
230

231
232
        ax.set_xlabel("mm")
        ax.set_ylabel("mm")
233
234
        # ax.set_title('edge_co')

235
        ax.legend(["edge_coords", "refined_edge_coords"])
236

237
        fig.savefig(f"{OUTDIR}/refined_coords.png", bbox_inches="tight", dpi=150)
238

239
    return refined_edge_coords
240
241


242
def _img_props(img, img_resolution, verbose=False):
243

244
245
    # # convert to grayscale
    # image = rgb2gray(img)
246

247
248
249
250
251
252
    # # find threshold for background
    # threshold_value = filters.threshold_otsu(image)

    image = rgb2hsv(img)[..., -1]
    image = exposure.equalize_hist(image)
    threshold_value = filters.threshold_li(image)
253
254

    if verbose:
255
        print(f"threshold={threshold_value}")
256
257
258
259
260
261
262
263
264
265
266
267
268
269

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

270
    bbox = np.array(p.bbox) * img_resolution
271
272

    centroid = (
273
274
        (bbox[0] + bbox[2]) / 2,
        (bbox[1] + bbox[3]) / 2,
275
276
    )

277
    if DO_PLOTS:
278

279
280
        fig, ax = plt.subplots()
        extent = np.array([0, img.shape[1], img.shape[0], 0]) * img_resolution
281

282
283
        im = ax.imshow(cc == p.label, extent=extent)
        fig.colorbar(im)
284

285
286
        rect = patches.Rectangle(
            (bbox[1], bbox[0]),
287
288
            bbox[3] - bbox[1],
            bbox[2] - bbox[0],
289
            linewidth=1,
290
291
            edgecolor="r",
            facecolor="none",
292
293
294
        )
        ax.add_patch(rect)

295
296
297
298
299
300
301
        ax.plot(centroid[1], centroid[0], "ro")

        ax.set_xlabel("mm")
        ax.set_ylabel("mm")
        ax.set_title("Image bounding box")

        fig.savefig(f"{OUTDIR}/image_bounding_box.png", bbox_inches="tight", dpi=150)
302

303
    return {
304
305
306
307
        "bbox": bbox,
        "bbox_centroid": centroid,
        "shape": (bbox[2] - bbox[0], bbox[3] - bbox[1]),
    }
308
309


310
def _pnt_props(pnts, verbose=False):
311
312
313
314

    x = pnts[:, 0]
    y = pnts[:, 1]

315
    bbox = (np.min(y), np.min(x), np.max(y), np.max(x))
316
317

    centroid = (
318
319
        (bbox[0] + bbox[2]) / 2,
        (bbox[1] + bbox[3]) / 2,
320
321
322
    )

    if verbose:
323
324
        print(f"bbox={bbox}")
        print(f"centroid={centroid}")
325

326
327
328
329
    if DO_PLOTS:

        fig, ax = plt.subplots()

330
331
        ax.plot(x, y, "b.")
        ax.axis("equal")
332
        ax.invert_yaxis()
333
334
335

        rect = patches.Rectangle(
            (bbox[1], bbox[0]),
336
337
            bbox[3] - bbox[1],
            bbox[2] - bbox[0],
338
            linewidth=1,
339
340
            edgecolor="r",
            facecolor="none",
341
        )
342
        ax.add_patch(rect)
343

344
345
346
        plt.plot(centroid[1], centroid[0], "ro")

        ax.set_title("Chart bounding box")
347

348
        fig.savefig(f"{OUTDIR}/chart_bounding_box.png", bbox_inches="tight", dpi=150)
349

350
351
352
353
354
    return {
        "bbox": bbox,
        "bbox_centroid": centroid,
        "shape": (bbox[2] - bbox[0], bbox[3] - bbox[1]),
    }
355
356


357
def _init_scale(img, img_resolution, crd, verbose=False):
358

359
    img_p = _img_props(img, img_resolution)
360
361
362
    crd_p = _pnt_props(crd)

    if verbose:
363
364
        print(f"image props = {img_p}")
        print(f"chart props = {crd_p}")
365
366
367

    idx = np.array([[1, 2], [1, 0], [3, 2], [3, 0]])

368
369
    src = np.array(crd_p["bbox"])[idx]
    dest = np.array(img_p["bbox"])[idx]
370

Sean Fitzgibbon's avatar
Sean Fitzgibbon committed
371
    a = transform.SimilarityTransform()
372
373
374
375
376
377
    a.estimate(src, dest)

    return a


def _apply_xfm(xfm, pnts):
378
379
380
    return (
        xfm.params @ np.concatenate((pnts, np.ones((pnts.shape[0], 1))), axis=-1).T
    ).T[:, :2]
381
382


383
def normal(pnts):
384
385
386
387
388
    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