Skip to content
Snippets Groups Projects
Commit ee51d512 authored by abda's avatar abda
Browse files

Solutions week 6

parent feb9ad49
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Week06/Plusplus_demo.ipynb 0 → 100644
+ 215
− 0
View file @ ee51d512
Source diff could not be displayed: it is too large. Options to address this: view the blob.
Week06/Plusplus_segmentation.ipynb 0 → 100644
+ 136
− 0
View file @ ee51d512
Source diff could not be displayed: it is too large. Options to address this: view the blob.
Week06/QUIZ_week6_solution.ipynb 0 → 100644
+ 131
− 0
View file @ ee51d512
%% Cell type:markdown id:26143b57 tags:
# QUIZ WEEK 6 SOLUTIONS
Questions:
* Inside area of the snake?
* Value of force affecting the point?
* Value of external energy?
%% Cell type:code id:97e7ba9e tags:
``` python
# Question 1
import numpy as np
print(f'inside_area: {np.pi * (100 + 10) ** 2}')
print(f'tolerance {(np.pi - 3.14) * (100 + 10) ** 2}')
```
%% Output
inside_area: 38013.2711084365
tolerance 19.271108436495197
%% Cell type:code id:d0c95345 tags:
``` python
# Questions 2 and 3
import numpy as np
import skimage
I = skimage.io.imread('../data/week6/plusplus.png', as_gray=True)
angles = np.linspace(0, 2*np.pi, 180, endpoint=False)
unit_circle = np.array([np.cos(angles), np.sin(angles)]).T
center = (np.array(I.shape).reshape([1, 2]) - 1) / 2
snake = center + 120 * unit_circle
mask = skimage.draw.polygon2mask(I.shape, snake)
m_in = (mask * I).sum()/mask.sum()
m_out =(~mask * I).sum()/(~mask).sum()
v = 0.8
f_ext = (m_in - m_out) * (2 * v-m_in - m_out)
print(f'{f_ext = }')
external_energy = (((I - m_in) ** 2) * mask).sum() + (((I - m_out) ** 2) * (~mask)).sum()
print(f'{external_energy = }')
```
%% Output
f_ext = 0.42836821376043044
external_energy = 23265.202043454898
%% Cell type:code id:82f41520 tags:
``` python
# Visualizing to make sure everything is correct
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.imshow(I, cmap=plt.cm.gray, vmin=0, vmax=1)
ax.imshow(mask, alpha=0.5)
ax.plot(snake[:, 0], snake[:, 1], 'r')
plt.show()
```
%% Output
Week06/simple_snake_new.py 0 → 100644
+ 122
− 0
View file @ ee51d512
#%%
"""
Representing snake as nx2 array of points.
@author: vand
"""
import numpy as np
import scipy.interpolate
import scipy.linalg
import skimage.draw
def make_circular_snake(N, center, radius):
""" Initialize circular snake."""
center = np.asarray(center).reshape([1, 2])
angles = np.linspace(0, 2*np.pi, N, endpoint=False)
unit_circle = np.array([np.cos(angles), np.sin(angles)]).T
return center + radius * unit_circle
def normalize(n):
l = np.sqrt((n ** 2).sum(axis=1, keepdims = True))
l[l == 0] = 1
return n / l
def get_normals(snake):
""" Returns snake normals. """
ds = normalize(np.roll(snake, 1, axis=0) - snake)
tangent = normalize(np.roll(ds, -1, axis=0) + ds)
normal = np.stack([-tangent[:, 1], tangent[:, 0]], axis=1)
return normal
def distribute_points(snake):
""" Distributes snake points equidistantly."""
N = len(snake)
closed = snake[np.hstack([np.arange(N), 0])]
d = np.sqrt(((np.roll(closed, 1, axis=0) - closed) ** 2).sum(axis=1))
d = np.cumsum(d)
d = d / d[-1] # Normalize to 0-1
x = np.linspace(0, 1, N, endpoint=False) # New points
new = np.stack([np.interp(x, d, closed[:, i]) for i in range(2)], axis=1)
return new
def is_ccw(A, B, C):
# Check if A, B, C are in counterclockwise order
return (C[1] - A[1]) * (B[0] - A[0]) > (B[1] - A[1]) * (C[0] - A[0])
def is_crossing(A, B, C, D):
# Check if line segments AB and CD intersect, not robust but ok for our case
# https://bryceboe.com/2006/10/23/line-segment-intersection-algorithm/
return is_ccw(A, C, D) != is_ccw(B, C, D) and is_ccw(A, B, C) != is_ccw(A, B, D)
def is_counterclockwise(snake):
""" Check if points are ordered counterclockwise."""
return np.dot(snake[1:, 0] - snake[:-1, 0],
snake[1:, 1] + snake[:-1, 1]) < 0
def remove_intersections(snake, method = 'new'):
""" Reorder snake points to remove self-intersections.
Arguments: snake represented by a N-by-2 array.
Returns: snake.
"""
N = len(snake)
closed = snake[np.hstack([np.arange(N), 0])]
for i in range(N - 2):
for j in range(i + 2, N):
if is_crossing(closed[i], closed[i + 1], closed[j], closed[j + 1]):
# Reverse vertices of smallest loop
rb, re = (i + 1, j) if j - i < N // 2 else (j + 1, i + N)
indices = np.arange(rb, re+1) % N
closed[indices] = closed[indices[::-1]]
snake = closed[:-1]
return snake if is_counterclockwise(snake) else np.flip(snake, axis=0)
def keep_snake_inside(snake, shape):
""" Contains snake inside the image."""
snake[:, 0] = np.clip(snake[:, 0], 0, shape[0] - 1)
snake[:, 1] = np.clip(snake[:, 1], 0, shape[1] - 1)
return snake
def regularization_matrix(N, alpha, beta):
""" Matrix for smoothing the snake."""
s = np.zeros(N)
s[[-2, -1, 0, 1, 2]] = (alpha * np.array([0, 1, -2, 1, 0]) +
beta * np.array([-1, 4, -6, 4, -1]))
S = scipy.linalg.circulant(s)
return scipy.linalg.inv(np.eye(N) - S)
def evolve_snake(snake, I, B, step_size):
""" Single step of snake evolution."""
mask = skimage.draw.polygon2mask(I.shape, snake)
m_in = np.mean(I[mask])
m_out = np.mean(I[~mask])
f = scipy.interpolate.RectBivariateSpline(np.arange(I.shape[0]), np.arange(I.shape[1]), I)
val = f(snake[:, 0], snake[:, 1], grid=False)
# val = I[snake[:, 0].astype(int), snake[:, 1].astype(int)] # simpler variant without interpolation
f_ext = (m_in - m_out) * (2 * val -m_in - m_out)
displacement = step_size * f_ext[:,None] * get_normals(snake)
snake = snake + displacement # external part
snake = B @ snake # internal part, ordering influenced by 2-by-N representation of snake
snake = remove_intersections(snake)
snake = distribute_points(snake)
snake = keep_snake_inside(snake, I.shape)
return snake
Week06/tracking_amoeba.ipynb 0 → 100644
+ 197
− 0
View file @ ee51d512
Source diff could not be displayed: it is too large. Options to address this: view the blob.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment