CellularPottsModel/cell_potts.py at master · lenkaback/CellularPottsModel · 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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
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
import numpy as np
import time
import cv2
import cmapy


# initialize the model and its parameters
class CPM:
    def __init__(self, beta=1., l=6, p=20000, q=1000, q1=40, tip_area=60, lambda_p=1, size=(600, 600), mu=0.02,
                 wound_area=15, div=1, fill=0):
        self.beta = beta  # inverse temperature
        self.l = l  # length of the sketch
        self.P = p  # number of spins tried to be updated per frame (update, proliferate, wound steps)
        self.q = q  # max amount of differently colored cells
        self.q1 = q1  # initial amount of differently colored cells
        self.counter = 0  # step counter of P updates
        self.J = np.empty((q + 1, q + 1))  # interaction strength matrix with similar color
        self.area = np.zeros(q)   # ideal area for each cell color
        self.tip_area = tip_area  # initial ideal area
        self.lambda_p = lambda_p  # area constraint strength
        self.colorz = np.zeros(q)  # current area of each cell (amount of pixels with the same color) - redundant?
        self.size = (self.height, self.width) = size  # if you change l, change here accordingly
        self.spin = np.zeros((self.width // self.l, self.height // self.l), dtype=np.uint8)  # lattice of cells
        self.mu = mu  # probability of proliferation
        self.wound_area = wound_area  # area of inflicted wound
        self.div = div  # should the cells proliferate? only at the beginning
        self.fill = fill  # is the lattice full of cells

        # initialization of all matrices
        # interactions: cell i/cell j (j_diff_cell), cell i/cell i (j_same_cell), cell/background (j_medium_cell)
        j_same_cell, j_diff_cell, j_medium_cell = 0, -5, -10
        self.J[:, :] = j_diff_cell
        for i in range(self.q + 1):
            self.J[0, i] = self.J[i, 0] = j_medium_cell
        np.fill_diagonal(self.J, j_same_cell)

        # all created cells have ideal area set at initial tip_area
        self.area[:q1] = self.tip_area

        # q1 initial cells created at random positions and updated current areas at 1 pixel
        for i in range(self.q1):
            i1 = np.random.randint(self.width // self.l)
            j1 = np.random.randint(self.height // self.l)

            self.spin[i1, j1] = i + 1
            self.colorz[i] = 1

    # try to make one spin update to change spin at position s(i2, j2) to s(i1, j1)
    def single_update(self):
        size_x, size_y = self.spin.shape
        # take random index of the updated position
        i1 = np.random.randint(size_x)
        j1 = np.random.randint(size_y)

        # pick a neighbor at random for position i2, j2
        ipi, jpi = np.random.choice([-1, 0, 1], 2)

        # take the neighboring position - continuous boundary condition
        i2 = (size_x + i1 + ipi) % size_x
        j2 = (size_y + j1 + jpi) % size_y

        # try to update the spin if not the same
        chosen = np.zeros(self.spin.shape)
        if self.spin[i1, j1] != self.spin[i2, j2]:
            # get all neighbors of (i1, j1) - 8-connectivity
            ip = (size_x + i1 + 1) % size_x
            im = (size_x + i1 - 1) % size_x
            jp = (size_y + j1 + 1) % size_y
            jm = (size_y + j1 - 1) % size_y

            square = 1 / np.sqrt(2)
            neighbors = np.array([[i1, jm], [i1, jp], [im, j1], [ip, j1]])
            neighbors_square = np.array([[im, jm], [im, jp], [ip, jm], [ip, jp]])

            chosen[:, :] = 0
            for px in neighbors:
                chosen[tuple(px)] = 1
            for px in neighbors_square:
                chosen[tuple(px)] = square

            neighbor_spins = self.spin.copy()
            neighbor_spins[chosen == 0] = 0

            spin1 = self.spin[i1, j1]
            spin2 = self.spin[i2, j2]

            F1 = np.sum([m * self.J[n, spin1] for n, m in zip(neighbor_spins[chosen > 0].reshape(-1), chosen[chosen > 0].reshape(-1))])
            F1n = np.sum([m * self.J[n, spin2] for n, m in zip(neighbor_spins[chosen > 0].reshape(-1), chosen[chosen > 0].reshape(-1))])


            # print("F1")
            # time_now = time.time()
            # F1 = self.J[spin1, self.spin[i1, jm]] + self.J[spin1, self.spin[i1, jp]]\
            #      + self.J[spin1, self.spin[im, j1]] + self.J[spin1, self.spin[ip, j1]]\
            #      + (self.J[spin1, self.spin[im, jm]] + self.J[spin1, self.spin[im, jp]]
            #      + self.J[spin1, self.spin[ip, jm]] + self.J[spin1, self.spin[ip, jp]]) / np.sqrt(2)
            # # I calculate the energy variation if I swap
            # time_f1 = time_now - time.time()
            # print(time_f1)
            # print("F1N")
            # time_now = time.time()
            # F1n = self.J[spin2, self.spin[i1, jm]] + self.J[spin2, self.spin[i1, jp]] \
            #     + self.J[spin2, self.spin[im, j1]] + self.J[spin2, self.spin[ip, j1]] \
            #     + (self.J[spin2, self.spin[im, jm]] + self.J[spin2, self.spin[im, jp]]
            #     + self.J[spin2, self.spin[ip, jm]] + self.J[spin2, self.spin[ip, jp]]) / np.sqrt(2)
            # time_f1n = time_f1 = time_now - time.time()
            # print(time_f1n)

            # print(F1, F1_my)
            # print(F1n, F1n_my)

            # calculate the difference in area term in old and new configuration, derivation of formulas in Overleaf
            deltaA = 0
            if self.spin[i1, j1] > 0 and self.spin[i2, j2] > 0:
                # deltaA = 2 * (np.sum(self.spin == self.spin[i2, j2]) - np.sum(self.spin == self.spin[i1, j1]) + 1)
                deltaA = 2 * (self.colorz[self.spin[i2, j2] - 1] - self.colorz[self.spin[i1, j1] - 1] + 1)
            else:
                if self.spin[i1, j1] == 0:
                    # deltaA = 1 + 2 * (np.sum(self.spin == self.spin[i2, j2]) - self.area[self.spin[i2, j2] - 1])
                    deltaA = 1 + 2 * (self.colorz[self.spin[i2, j2] - 1] - self.area[self.spin[i2, j2] - 1])
                else:
                    # deltaA = 1 - 2 * (np.sum(self.spin == self.spin[i1, j1]) - self.area[self.spin[i1, j1] - 1])
                    deltaA = 1 - 2 * (self.colorz[self.spin[i1, j1] - 1] - self.area[self.spin[i1, j1] - 1])
            delta = F1 - F1n + self.lambda_p * deltaA  # energy variation of the swap

            # Metropolis rule of swap (if lower energy or if Boltzmann condition is satisfied)
            cas = np.random.rand(1)
            if delta <= 0 or cas <= 1 / np.exp(self.beta * delta):
                news = 1  # Metropolis rule
                if self.spin[i1, j1] > 0:
                    self.colorz[self.spin[i1, j1]-1] -= 1
                if self.spin[i2, j2] > 0:
                    self.colorz[self.spin[i2, j2]-1] += 1
                self.spin[i1, j1] = self.spin[i2, j2]

    def initial_proliferation(self):
        cas2 = np.random.rand(1)
        if cas2 < self.mu:
            size_x, size_y = self.spin.shape
            i3 = np.random.randint(size_x)
            j3 = np.random.randint(size_y)
            if self.spin[i3, j3] > 0:
                self.area[self.spin[i3, j3] - 1] += 1
                if self.area[self.spin[i3, j3] - 1] == 2 * self.tip_area and self.div == 1:
                    self.area[self.spin[i3, j3] - 1] = self.tip_area
                    q2 = int(np.random.rand(1) * self.q)
                    while self.area[q2] > 0:
                        q2 = int(np.random.rand(1) * self.q)

                    self.spin[i3, j3] = q2 + 1
                    self.area[q2] = self.tip_area
                    self.colorz[q2] += 1
                    self.colorz[self.spin[i3, j3] - 1] -= 1

    def wound_formation(self):
        if self.fill == 1 and self.counter % 200 == 0:
            for u in range(self.q):
                if self.area[u] > 0:
                    print("Wounding")
                    cent_massx = 0
                    cent_massy = 0
                    for i4 in range(self.width // self.l):
                        for j4 in range(self.height // self.l):
                            if self.spin[i4, j4] == u + 1:
                                cent_massx += i4
                                cent_massy += j4

                    cent_massx /= self.colorz[u]
                    cent_massy /= self.colorz[u]
                    dist = np.sqrt((cent_massx - 50) * (cent_massx - 50)
                                + (cent_massy - 50) * (cent_massy - 50))
                    if dist < self.wound_area:
                        self.area[u] = self.colorz[u] = 0
                        for i4 in range(self.width // self.l):
                            for j4 in range(self.height // self.l):
                                if self.spin[i4, j4] == u+1:
                                    self.spin[i4, j4] = 0


np.random.seed(42)
cpm = CPM()
cv2.namedWindow('img', cv2.WINDOW_NORMAL)
# exit()
print("all")
time_now = time.time()
for counter in range(5000):
    for i in range(cpm.P):
        # print("update")
        # time_now = time.time()
        cpm.single_update()
        # time_upd = time.time() - time_now
        # print(time_upd)
        #
        # print("prolif")
        # time_now = time.time()
        cpm.initial_proliferation()
        # time_upd = time.time() - time_now
        # print(time_upd)
        cpm.wound_formation()

        if np.sum(cpm.spin > 0) >= 0.99 * cpm.width * cpm.height / (cpm.l * cpm.l):
            if cpm.fill != 1:
                print('FULL')
                cpm.fill = 1
                cpm.mu = 0
                cpm.div = 0
                cpm.lambda_p = 0.2
        else:
            cpm.fill = 0
    cpm.counter += 1

    lattice_cm = cv2.applyColorMap(12 * cpm.spin.astype(np.uint8), cmapy.cmap('gist_ncar_r'))
    cv2.imshow('img', lattice_cm)
    cv2.waitKey(1)
time_upd = time.time() - time_now
print(time_upd)