Skip to content
Snippets Groups Projects
Commit 6cc80b24 authored by Vedrana Andersen Dahl's avatar Vedrana Andersen Dahl
Browse files

Added week 5 solutions

parent 309c1196
No related branches found
No related tags found
No related merge requests found
%% Cell type:markdown id: tags:
# Binary segmentation of DTU image
%% Cell type:code id: tags:
``` python
from skimage.io import imread
import matplotlib.pyplot as plt
import maxflow
import numpy as np
def segmentation_histogram(ax, I, S, edges=None):
'''
Histogram for data and each segmentation label.
'''
if edges is None:
edges = np.linspace(I.min(), I.max(), 100)
ax.hist(I.ravel(), bins=edges, color = 'k')
centers = 0.5*(edges[:-1] + edges[1:]);
for k in range(S.max()+1):
ax.plot(centers, np.histogram(I[S==k].ravel(), edges)[0])
```
%% Cell type:code id: tags:
``` python
# noisy image
I = imread('../data/week5/DTU_noisy.png').astype(float)/255
# MRF parameters
beta = 0.1
mu = [90/255, 170/255]
# Setting up graph with internal and external edges
g = maxflow.Graph[float]()
nodeids = g.add_grid_nodes(I.shape)
g.add_grid_edges(nodeids, beta)
g.add_grid_tedges(nodeids, (I - mu[1])**2, (I - mu[0])**2)
# Graph cut
g.maxflow()
S = g.get_grid_segments(nodeids)
# Visualization
fig, ax = plt.subplots(1, 3)
ax[0].imshow(I, vmin=0, vmax=1, cmap=plt.cm.gray)
ax[0].set_title('Noisy image')
ax[1].imshow(S)
ax[1].set_title('Segmented')
segmentation_histogram(ax[2], I, S, edges=None)
ax[2].set_aspect(1./ax[2].get_data_ratio())
ax[2].set_title('Segmentation histogram')
plt.tight_layout()
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id:26143b57 tags:
# QUIZ WEEK 5 SOLUTIONS
%% Cell type:markdown id:63a90f84 tags:
## Questions:
- What is the prior energy of the maximum likelihood solution?
- What is the likelihood energy of the row-configuration?
- What is the posterior energy of column-configuration?
%% Cell type:code id:d0c95345 tags:
``` python
import numpy as np
def prior_energy(S, beta):
return beta * ((S[:-1] != S[1:]).sum() + (S[:, :-1] != S[:, 1:]).sum())
def likelihood_energy(S, I, mu):
return ((mu[S] - I)**2).sum()
I = np.array([[1, 2, 6, 4, 10, 8],
[4, 1, 3, 5, 9, 7],
[5, 2, 3, 5, 4, 8]])
mu = np.array([2, 5, 10])
beta = 10
#%% Question 1
U = (np.stack([I - mu[0], I - mu[1], I - mu[2]], axis=2) ** 2)
S0 = np.argmin(U, axis=2) # max likelihood solution
prior_max_likelihood = prior_energy(S0, beta)
print(f'{prior_max_likelihood = }')
#%% Question 2
S = np.tile([[0], [1], [2]], (1, 6)).astype(int)
likelihood_rows = likelihood_energy(S, I, mu)
print(f'{likelihood_rows = }')
#%% Question 3
S = np.tile([0, 0, 1, 1, 2, 2], (3, 1)).astype(int)
prior_stripes = prior_energy(S, beta)
likelihood_stripes = likelihood_energy(S, I, mu)
posterior_stripes = prior_stripes + likelihood_stripes
print(f'{posterior_stripes = }')
#%% Extra solution, not asked in the quiz - computing MAP using graph cuts
import maxflow.fastmin
S_GC = S0.copy()
maxflow.fastmin.aexpansion_grid(U, beta -
beta*np.eye(3, 3, dtype=U.dtype), labels = S_GC)
posterior_MAP = prior_energy(S_GC, beta) + likelihood_energy(S_GC, I, mu)
print(f'EXTRA, {posterior_MAP = }')
# %%
```
%% Output
prior_max_likelihood = 140
likelihood_rows = 365
posterior_stripes = 139
EXTRA, posterior_MAP = 114
%% Cell type:code id:cc2643c4 tags:
``` python
```
This diff is collapsed.
Source diff could not be displayed: it is too large. Options to address this: view the blob.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment