2dShapeSpacePortable/shapespace.py at main · CellProfiling/2dShapeSpacePortable · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
import numpy as np
from scipy.ndimage import center_of_mass
from skimage.measure import find_contours
import matplotlib.pyplot as plt
import os

from lib import helpers, alignment, image_utils, coefs


def shapespace(config, image_id, nuc_path, cell_maskpath, protein_path):
    nuc_mask = image_utils.read_grayscale_image(nuc_path)
    cell_mask = image_utils.read_grayscale_image(cell_maskpath)
    prot_img = image_utils.read_grayscale_image(protein_path)

    nuc_area = np.sum(nuc_mask > 0)
    cell_area = np.sum(cell_mask > 0)
    prot_nuc_sum = np.sum(prot_img * nuc_mask)
    prot_cell_sum = np.sum(prot_img * cell_mask)

    if config['alignment'] == "fft_major_axis":
        nuc_mask, cell_mask, theta = alignment.align_major_axis(nuc_mask, cell_mask, prot_img, plot=False)
    elif config['alignment'] == "fft_major_axis_polarized":
        nuc_mask, cell_mask, theta = alignment.align_major_axis_polarized(nuc_mask, cell_mask, prot_img, plot=False)
    elif config['alignment'] == "fft_centroid":
        nuc_mask, cell_mask, theta = alignment.align_centroids(nuc_mask, cell_mask, prot_img, plot=False)
    else:
        raise NotImplementedError(f"Unknown alignment method: {config['alignment']}")
    centroid = center_of_mass(nuc_mask)

    # Padd surrounding with 0 so no contour touch the border. This help matching squares algo not failing (as much)
    nuclei = np.zeros((nuc_mask.shape[0] + 2, nuc_mask.shape[1] + 2))
    nuclei[1 : 1 + nuc_mask.shape[0], 1 : 1 + nuc_mask.shape[1]] = nuc_mask
    cell = np.zeros((cell_mask.shape[0] + 2, cell_mask.shape[1] + 2))
    cell[1 : 1 + cell_mask.shape[0], 1 : 1 + cell_mask.shape[1]] = cell_mask

    nuc_maskcoords_ = find_contours(nuclei)
    if (
        len(nuc_maskcoords_) > 1
    ):  # concatenate fragmented contour lines, original point could be ambiguous! (attempt to re-align original point in coefs.XXX_fourier_coefs())
        # if the biggest contour segment is a close loop, the rest are micronuclei, artifact segments
        idx_longest = np.argmax([len(xy) for xy in nuc_maskcoords_])
        biggest = nuc_maskcoords_[idx_longest]
        if all(biggest[0] == biggest[-1]):
            nuc_maskcoords_ = biggest - centroid
        else:
            nuc_maskcoords_ = np.vstack(nuc_maskcoords_)
            nuc_maskcoords_ = nuc_maskcoords_ - centroid
    else:
        nuc_maskcoords_ = nuc_maskcoords_[0] - centroid

    cell_maskcoords_ = find_contours(cell)
    if (
        len(cell_maskcoords_) > 1
    ):  # concatenate fragmented contour lines, original point could be ambiguous! (attempt to re-align original point in coefs.XXX_fourier_coefs())
        # if the biggest contour segment is a close loop, the rest are micronuclei, artifact segments
        idx_longest = np.argmax([len(xy) for xy in cell_maskcoords_])
        biggest = cell_maskcoords_[idx_longest]
        if all(biggest[0] == biggest[-1]):
            cell_maskcoords_ = biggest - centroid
        else:
            cell_maskcoords_ = np.vstack(cell_maskcoords_)
            cell_maskcoords_ = cell_maskcoords_ - centroid
    else:
        cell_maskcoords_ = cell_maskcoords_[0] - centroid

    cell_maskcoords = cell_maskcoords_.copy()
    cell_maskcoords = helpers.realign_contour_startpoint(cell_maskcoords)
    nuc_maskcoords = nuc_maskcoords_.copy()
    nuc_maskcoords = helpers.realign_contour_startpoint(nuc_maskcoords)
    if config["plot"]:
        fig, ax = plt.subplots(1, 3, figsize=(8, 4))
        ax[0].imshow(nuc_mask, alpha=0.5, origin="lower")
        ax[0].imshow(cell_mask, alpha=0.5, origin="lower")

        nu_centroid = helpers.find_centroid(nuc_maskcoords)
        cell_maskcentroid = helpers.find_centroid(cell_maskcoords)
        ax[1].plot(nuc_maskcoords_[:, 0], nuc_maskcoords_[:, 1])
        ax[1].scatter(nuc_maskcoords_[0, 0], nuc_maskcoords_[0, 1], color="slateblue")
        ax[1].scatter(nu_centroid[0], nu_centroid[1], color="b")
        ax[1].plot(cell_maskcoords_[:, 0], cell_maskcoords_[:, 1])
        ax[1].scatter(cell_maskcoords_[0, 0], cell_maskcoords_[0, 1], color="gold")
        ax[1].scatter(cell_maskcentroid[0], cell_maskcentroid[1], color="orange")

        ax[1].vlines(0, -200, 200, colors="gray", linestyles="dashed")
        ax[1].hlines(0, -200, 200, colors="gray", linestyles="dashed")
        ax[1].axis("scaled")
        ax[2].set_title(f"theta = {np.round(theta,1)}°")
        ax[2].vlines(0, -200, 200, colors="gray", linestyles="dashed")
        ax[2].hlines(0, -200, 200, colors="gray", linestyles="dashed")
        ax[2].plot(nuc_maskcoords[:, 0], nuc_maskcoords[:, 1])
        ax[2].scatter(nuc_maskcoords[0, 0], nuc_maskcoords[0, 1], color="slateblue")
        ax[2].scatter(nu_centroid[0], nu_centroid[1], color="b")
        ax[2].plot(cell_maskcoords[:, 0], cell_maskcoords[:, 1])
        ax[2].scatter(cell_maskcoords[0, 0], cell_maskcoords[0, 1], color="gold")
        ax[2].scatter(cell_maskcentroid[0], cell_maskcentroid[1], color="orange")
        ax[2].axis("scaled")
        plt.savefig(os.path.join(config["output_dir"], "shapespace", image_id + "_alignment.png"), bbox_inches="tight")
        plt.close()

    fcoef_n, e_n = coefs.fourier_coeffs(nuc_maskcoords, n_coef=config["n_coeffs"], plot=config["plot"],
                                         save_path=os.path.join(config["output_dir"], "shapespace", image_id + "_fft_nuc.png"))
    fcoef_c, e_c = coefs.fourier_coeffs(cell_maskcoords, n_coef=config["n_coeffs"], plot=config["plot"],
                                         save_path=os.path.join(config["output_dir"], "shapespace", image_id + "_fft_cell.png"))
    final_list = [image_id, nuc_area, cell_area, prot_nuc_sum, prot_cell_sum, theta, str(centroid[0]), str(centroid[1]), e_c, e_n]
    final_list.extend(fcoef_n[0])
    final_list.extend(fcoef_n[1])
    final_list.extend(fcoef_c[0])
    final_list.extend(fcoef_c[1])
    return [str(element) for element in final_list]