Skip to content
Snippets Groups Projects
Commit c3f3971b authored by RTuxen's avatar RTuxen
Browse files

Added describing text

parent 5c146221
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Showing
with 1218 additions and 36 deletions
notebooks/.ipynb_checkpoints/LayerSegmentation2D-checkpoint.ipynb 0 → 100644
+ 321
− 0
View file @ c3f3971b
%% Cell type:markdown id: tags:
# 2D Layer Surface Segmentation
%% Cell type:markdown id: tags:
**Add Proper References Here**
%% Cell type:markdown id: tags:
Load essential modules for loading and showing data
%% Cell type:code id: tags:
``` python
# Load modules
import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
```
%% Cell type:code id: tags:
``` python
# Load image
in_dir = 'data/'
data = imread(in_dir+'layer2D.png').astype(np.int32)
# Show image.
plt.imshow(data, cmap='gray')
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
## Detect one layer
%% Cell type:code id: tags:
``` python
from slgbuilder import GraphObject,MaxflowBuilder
```
%% Cell type:markdown id: tags:
Let's create an object that represent the layer we're trying a segment. Here we are only segmenting one layer. The image data is passed directly to the ```GraphObject```. By default, it assumes data is structured on a regular index grid.
The algorithm segments layers by detecting dark lines. For better segmentation the image is inversed by `255-data`.
%% Cell type:code id: tags:
``` python
layer_1 = GraphObject(255-data)
```
%% Cell type:markdown id: tags:
Next we create the ```MaxflowBuilder``` and add the layer. Then we add layered cost and smoothness. Edges are added automatically for when calling ```add_layered_boundary_cost``` and ```add_layered_smoothness```. By default these edges are added for all object previously added to the ```GraphHelper```. However, it is possible to add the edges for specific object only, by passing a list of objects to the function. The graph structure is based on [Li et al](https://doi.org/10.1109/TPAMI.2006.19).
%% Cell type:code id: tags:
``` python
helper = MaxflowBuilder()
helper.add_object(layer_1)
helper.add_layered_boundary_cost()
helper.add_layered_smoothness()
```
%% Cell type:markdown id: tags:
The graph is cut using the [Maxflow](https://github.com/Skielex/maxflow) algorithm to get the maximum flow/minimum cut. This algorithm will find the global optimal solution to the binary segmentation problem.
%% Cell type:code id: tags:
``` python
flow = helper.solve()
print('Maximum flow/minimum energy:', flow)
```
%% Output
Maximum flow/minimum energy: 86761
%% Cell type:markdown id: tags:
The segmentation for ```layer_1``` can then be displayed.
%% Cell type:code id: tags:
``` python
segmentation = helper.what_segments(layer_1)
segmentation_line = np.argmin(segmentation, axis=0)
# Draw results.
f,ax = plt.subplots(1,3,figsize=(10,10))
ax[0].imshow(data, cmap='gray')
ax[1].imshow(segmentation)
ax[2].imshow(data, cmap='gray')
ax[2].plot(segmentation_line)
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
## Detect multiple layers
%% Cell type:markdown id: tags:
This time we will create two objects, as we want to segment two different layers. Although our graph cut is still binary, we will create the graph so that we get a multi-label segmentation.
%% Cell type:code id: tags:
``` python
layers = []
for i in range(2):
layers.append(GraphObject(255-data))
```
%% Cell type:markdown id: tags:
We create the ```MaxflowBuilder``` and add all layers to the graph. We then add boundary costs and smoothness as before.
%% Cell type:code id: tags:
``` python
helper = MaxflowBuilder()
helper.add_objects(layers)
helper.add_layered_boundary_cost()
helper.add_layered_smoothness()
```
%% Cell type:markdown id: tags:
We've added the layers to the graph, but because they all are based on the same data, and there are no interaction contraints between the layers, they would all result in the exact same segmentation. This is *not* what we want.
We want the layers to be separated by a certain distance. Luckily [Li et al](https://doi.org/10.1109/TPAMI.2006.19). describe how to enforce a minimum distance between the layers. This technique has been implemented in ```add_layered_containment```, which takes an *outer* and an *inner* ```GraphObject``` and lets you specify both a minimum and maximum distance/margin between them. In this case we will only specify a minimum margin.
%% Cell type:code id: tags:
``` python
for i in range(len(layers)-1):
helper.add_layered_containment(layers[i], layers[i+1], min_margin=10)
```
%% Cell type:code id: tags:
``` python
flow = helper.solve()
print('Maximum flow/minimum energy:', flow)
```
%% Output
Maximum flow/minimum energy: 174018
%% Cell type:markdown id: tags:
We can now get the segmentations for all layers and show them. Notice that the blue line, representing the first layers is the furthest down, or outermost.
%% Cell type:code id: tags:
``` python
segmentations = [helper.what_segments(l).astype(np.int32) for l in layers]
segmentation_lines = [np.argmin(s, axis=0) - 0.5 for s in segmentations]
# Draw results.
f,ax = plt.subplots(1,3,figsize=(10,10))
ax[0].imshow(data, cmap='gray')
ax[1].imshow(np.sum(segmentations, axis=0))
ax[2].imshow(data, cmap='gray')
for line in segmentation_lines:
ax[2].plot(line)
plt.show()
```
%% Output
notebooks/.ipynb_checkpoints/LayerSegmentation3D-checkpoint.ipynb 0 → 100644
+ 207
− 0
View file @ c3f3971b
This diff is collapsed.
notebooks/LayerSegmentation2D.ipynb 0 → 100644
+ 416
− 0
View file @ c3f3971b
Source diff could not be displayed: it is too large. Options to address this: view the blob.
notebooks/LayerSegmentation3D.ipynb 0 → 100644
+ 237
− 0
View file @ c3f3971b
This diff is collapsed.
notebooks/NerveSegmentation2D.ipynb
+ 12
− 12
View file @ c3f3971b
%% Cell type:markdown id: tags:
 
# 2D Nerve Surface Segmentation
 
%% Cell type:markdown id: tags:
 
**Add Proper References Here**
 
%% Cell type:code id: tags:
 
``` python
# Load modules
import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
```
 
%% Cell type:markdown id: tags:
 
The objects (nerves) in the image have two different sizes which need to be handled differently. With small nerves marked in red an large in blue.
The objects (nerves) in the image have two different sizes which need to be handled differently. With small nerves marked in red and large nerves marked in blue.
 
%% Cell type:code id: tags:
 
``` python
# Load data
in_dir = 'data/'
data = imread(in_dir+'nerves2D.png').astype(np.int32)
 
# Get centers
path = in_dir+'nerveCenters.png'
centers = imread(path)[...,0]
 
centers_small = np.transpose(np.where(centers==1))
centers_large = np.transpose(np.where(centers==2))
 
 
# Show image with centers.
plt.imshow(data, cmap='gray')
plt.scatter(centers_small[..., 1], centers_small[..., 0], color='red', s=6)
plt.scatter(centers_large[..., 1], centers_large[..., 0], color='blue', s=6)
plt.show()
 
print(f'Number of objects: {len(centers_small)+len(centers_large)}')
```
 
%% Output
 
 
Number of objects: 25
 
%% Cell type:markdown id: tags:
 
## Unfolding
To detect the objects (nerves) using layered surface detection, we first need to unfold the nerves using a radial resampling.
 
%% Cell type:code id: tags:
 
``` python
from scipy.ndimage.interpolation import map_coordinates
 
def unfold_image(img, center, max_dists=None, r_min=1, r_max=20, angles=30, steps=15):
""" Unfolds around a point in an image with radial resampling """
# Sampling angles and radii.
angles = np.linspace(0, 2*np.pi, angles, endpoint=False)
distances = np.linspace(r_min, r_max, steps, endpoint=True)
 
if max_dists is not None:
max_dists.append(np.max(distances))
 
# Get angles.
angles_cos = np.cos(angles)
angles_sin = np.sin(angles)
 
# Calculate points positions.
x_pos = center[0] + np.outer(angles_cos, distances)
y_pos = center[1] + np.outer(angles_sin, distances)
 
# Create list of sampling points.
sampling_points = np.array([x_pos, y_pos]).transpose()
sampling_shape = sampling_points.shape
sampling_points_flat = sampling_points.reshape((-1, 2))
 
# Sample from image.
samples = map_coordinates(img, sampling_points_flat.transpose(), mode='nearest')
samples = samples.reshape(sampling_shape[:2])
 
return samples, sampling_points
```
 
%% Cell type:markdown id: tags:
 
Now that we have a function for unfolding image data, let's test it. The result should be an unfolded image, for which we can use layer detection.
 
%% Cell type:code id: tags:
 
``` python
samples, sample_points = unfold_image(data, centers_small[4],r_max=35,angles=30,steps=30)
 
plt.figure(figsize=(13, 5))
ax = plt.subplot(1, 3, 1, title='Sample positions in data')
ax.imshow(data, cmap='gray')
ax.scatter(sample_points[..., 1], sample_points[..., 0], s=2, color='red')
ax = plt.subplot(1, 3, 2, title='Sample positions and intensities')
ax.scatter(sample_points[..., 1], sample_points[..., 0], c=samples, cmap='gray')
ax = plt.subplot(1, 3, 3, title='Unfolded image')
ax.imshow(samples, cmap='gray')
plt.show()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
## Detect layers in object
Now that we can unfold the nerves, we can try to use graph cut based layer detection.
 
Since we want to separate the inner and outer part of the nerve, we will detect two layers per nerve. We will use the gradient image for this.
 
%% Cell type:code id: tags:
 
``` python
from slgbuilder import GraphObject,MaxflowBuilder
 
# Create gradient-based objects.
diff_samples = np.diff(samples, axis=0)
outer_nerve = GraphObject(255 - diff_samples)
inner_nerve = GraphObject(diff_samples)
 
f,ax = plt.subplots(1,2,figsize=(10,5))
ax[0].imshow(outer_nerve.data, cmap='gray')
ax[0].set_title('Outer never data')
ax[1].imshow(inner_nerve.data, cmap='gray')
ax[1].set_title('Inner never data')
plt.show()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
The surface will be detected where the data pixel intensity in the images are low. This corresponds well with the data for the outer and inner nerves shown above.
 
Let's detect the layers. We apply boundary cost, smoothness and containment constraints. Here we set both ```min_margin``` and ```max_margin``` constraints for our containment. Then we use ```maxflow``` to find the optimal solution.
 
%% Cell type:code id: tags:
 
``` python
helper = MaxflowBuilder()
helper.add_objects([outer_nerve, inner_nerve])
helper.add_layered_boundary_cost()
helper.add_layered_smoothness(delta=2)
helper.add_layered_containment(outer_nerve, inner_nerve, min_margin=3, max_margin=6)
 
flow = helper.solve()
print('Maximum flow/minimum energy:', flow)
```
 
%% Output
 
Maximum flow/minimum energy: 4754
 
%% Cell type:markdown id: tags:
 
We can then get the segmentations and the segmentation lines from the *GraphObject*
 
%% Cell type:code id: tags:
 
``` python
segmentations = [helper.get_labels(o).astype(np.int32) for o in helper.objects]
segmentation_lines = [np.count_nonzero(s, axis=0) - 0.5 for s in segmentations]
 
# Draw segmentation
f,ax = plt.subplots(1,3,figsize = (10,10))
ax[0].imshow(samples, cmap='gray')
ax[1].imshow(np.sum(segmentations,axis=0))
ax[2].imshow(samples, cmap='gray')
for line in segmentation_lines:
ax[2].plot(line)
plt.show()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
## Detecting multiple objects
In the image data here we have marked 25 different nerves that we would like to segment. We could segment each of these individually, the same way we segmented the single nerve above. To get a better segmentation the large nerves will need to be samples with slightly different parameters than the smaller ones. Although it is not the most memory efficient way of segmenting the objects, we could just add all the object to the graph at once and get a segmentation for each object.
 
%% Cell type:code id: tags:
 
``` python
# Lists for storing nerve objects.
nerve_samples = []
outer_nerves = []
inner_nerves = []
 
# For each center, create an inner and outer never.
for center in centers_small:
samples, sample_points = unfold_image(data, center,r_max=35,angles=40,steps=30)
nerve_samples.append(samples)
 
# Create outer and inner nerve objects.
diff_samples = np.diff(samples, axis=0)
diff_sample_points = sample_points[:-1]
 
 
outer_nerves.append(GraphObject(255 - diff_samples, diff_sample_points))
inner_nerves.append(GraphObject(diff_samples, diff_sample_points))
for center in centers_large:
samples, sample_points = unfold_image(data, center,r_max=60,angles=40,steps=30)
nerve_samples.append(samples)
 
# Create outer and inner nerve objects.
diff_samples = np.diff(samples, axis=0)
diff_sample_points = sample_points[:-1]
 
 
outer_nerves.append(GraphObject(255 - diff_samples, diff_sample_points))
inner_nerves.append(GraphObject(diff_samples, diff_sample_points))
```
 
%% Cell type:markdown id: tags:
 
Here we also add the sample positions to the ```GraphObject```s as they are needed to draw the sementations on the nerve image.
 
%% Cell type:code id: tags:
 
``` python
helper = MaxflowBuilder()
helper.add_objects(outer_nerves + inner_nerves)
helper.add_layered_boundary_cost()
helper.add_layered_smoothness(delta=2)
 
for outer_nerve, inner_nerve in zip(outer_nerves, inner_nerves):
helper.add_layered_containment(outer_nerve, inner_nerve, min_margin=3, max_margin=6)
 
flow = helper.solve()
print('Maximum flow/minimum energy:', flow)
```
 
%% Output
 
Maximum flow/minimum energy: 212430
 
%% Cell type:code id: tags:
 
``` python
# Get segmentations.
segmentations = []
for outer_nerve, inner_nerve in zip(outer_nerves, inner_nerves):
segmentations.append(helper.get_labels(outer_nerve))
segmentations.append(helper.get_labels(inner_nerve))
 
segmentation_lines = [np.count_nonzero(s, axis=0) - 0.5 for s in segmentations]
 
# Draw segmentations.
plt.figure(figsize=(15, 5))
for i, samples in enumerate(nerve_samples):
ax = plt.subplot(3, len(nerve_samples) // 3 + 1, i + 1)
ax.imshow(samples, cmap='gray')
 
ax.plot(segmentation_lines[2*i])
ax.plot(segmentation_lines[2*i + 1])
 
plt.show()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
While most of the segmentations went well, if we look closely we see that some don't look right. If we draw the lines on the original image, we see the problem.
 
%% Cell type:code id: tags:
 
``` python
def draw_segmentations(data, helper):
"""Draw all segmentations for objects in the helper on top of the data."""
 
# Create figure.
plt.figure(figsize=(10, 10))
plt.imshow(data, cmap='gray')
plt.xlim([0, data.shape[1]-1])
plt.ylim([data.shape[0]-1, 0])
 
# Draw segmentation lines.
for i, obj in enumerate(helper.objects):
 
# Get segmentation.
segment = helper.get_labels(obj)
 
# Create line.
line = np.count_nonzero(segment, axis=0)
 
# Get actual points.
point_indices = tuple(np.asarray([line - 1, np.arange(len(line))]))
points = obj.sample_points[point_indices]
# Close line.
points = np.append(points, points[:1], axis=0)
 
# Plot points.
plt.plot(points[..., 1], points[..., 0])
 
plt.show()
 
draw_segmentations(data, helper)
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
<!-- Some of the objects' segmentation are overlapping -->
 
%% Cell type:markdown id: tags:
 
Some of the objects' segmentation are overlapping
 
%% Cell type:markdown id: tags:
 
## Multi-object exclusion
To overcome the issue of overlapping segments, we can add exclusion contraints between all outer nerves. However, exclusion is a so-called *nonsubmodular* energy term, which means it cannot be represented as a single edge in our graph. Luckily there's an algorithm called *QPBO* (Qudratic Pseudo-Boolean Optimization) that can help us.
 
QPBO creates creates a complementary graph, alongside the original graph. The complementary graph is inverted, meaning that is has the exact same edges as the original graph, except they are reversed. This means that the graph size is doubled, which makes computation slower and uses more memory. The benefit of QPBO is that we can now add nonsubmodular energies such as exclusion.
 
The ```graphhelper``` module contains a ```QPBOHelper``` class, which is very similar to the ```GraphHelper``` we've been using so far. The main difference is that it has functions for adding exclusion. One of these is ```add_layered_exclusion``` which we will now use. We will be using the ```GraphObject```s created earlier.
 
%% Cell type:code id: tags:
 
``` python
from slgbuilder import QPBOBuilder
 
helper = QPBOBuilder()
helper.add_objects(outer_nerves + inner_nerves)
helper.add_layered_boundary_cost()
helper.add_layered_smoothness(delta=2)
 
for outer_nerve, inner_nerve in zip(outer_nerves, inner_nerves):
helper.add_layered_containment(outer_nerve, inner_nerve, min_margin=3, max_margin=6)
 
twice_flow = helper.solve()
print('Two times maximum flow/minimum energy:', twice_flow)
 
if 2*flow == twice_flow:
print('QPBO flow is exactly twice the Maxflow flow.')
else:
print('Something is wrong...')
```
 
%% Output
 
Two times maximum flow/minimum energy: 424860
QPBO flow is exactly twice the Maxflow flow.
 
%% Cell type:markdown id: tags:
 
We see that the ```QPBOHelper``` energy/flow is exactly twice the flow computed by ```GraphHelper``` for a similar problem, which is what we expect, since we double the number of nodes and edges. This is because we have added exactly the same edges/energies on above. This of course also means that the segmentation is exactly the same, hence we haven't fixed the problem yet.
 
To avoid the overlapping nerve segments, we add exclusion between all *outer* nerve objects using ```add_layered_exclusion``` and call ```solve``` again. Note that calculating the new maxflow/mincut only requires us to re-evaluate parts of the graph that were changed, potentially making the computation very fast.
 
%% Cell type:code id: tags:
 
``` python
# Add exclusion constraints between all pairs of outer nerves.
for i in range(len(outer_nerves)):
for j in range(i + 1, len(outer_nerves)):
helper.add_layered_exclusion(outer_nerves[i], outer_nerves[j], margin=3)
 
twice_flow = helper.solve()
print('Two times maximum flow/minimum energy:', twice_flow)
```
 
%% Output
 
Two times maximum flow/minimum energy: 425476
 
%% Cell type:markdown id: tags:
 
We see that adding the new constraints has increased the energy. This makes sense, since our constraints are forcing a solution that is less optimal from the perspective of the data. However, our prior knowledge tells us that nerves cannot overlap, so even if the data suggest that they do, we know this is not the case, but rather because the data is inaccurate.
 
Let's draw the segmentation results with exclusion inteactions.
 
%% Cell type:code id: tags:
 
``` python
draw_segmentations(data, helper)
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
Which yields a better overall segmentation
......
notebooks/NerveSegmentation3D.ipynb
+ 25
− 24
View file @ c3f3971b
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment