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

Ignore .py help-files

parent 69fa9f15
No related branches found
No related tags found
No related merge requests found
"""
This script does the exact same as 'NerveSegmenetation2D.ipynb'
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
from scipy.ndimage.interpolation import map_coordinates
from slgbuilder import GraphObject,MaxflowBuilder
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()
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 sampling """
# 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
#%%
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)}')
#%% Example of radial unfolding
samples, sample_points = unfold_image(data, centers_small[4],r_max=20,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()
#%% Detect layers in object
# 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()
#%% Segment one sample
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)
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]
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()
#%% Multiple object
# 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))
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)
# 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()
draw_segmentations(data, helper)
#%% Multi-object exclusion
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...')
# 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)
draw_segmentations(data, helper)
#%% Region cost
# mu_inside = 90
# mu_ring = 70
# mu_outside = 90
# beta = 0.1
# for samples, outer_nerve, inner_nerve in zip(nerve_samples, outer_nerves, inner_nerves):
# samples = samples[:-1]
# inside_cost = np.abs(samples - mu_inside) * beta
# ring_cost = np.abs(samples - mu_ring) * beta
# outside_cost = np.abs(samples - mu_outside) * beta
# helper.add_layered_region_cost(inner_nerve, ring_cost, inside_cost)
# helper.add_layered_region_cost(outer_nerve, outside_cost, ring_cost)
# twice_flow = helper.solve()
# print('Two times maximum flow/minimum energy:', twice_flow)
# draw_segmentations(data, helper)
"""
This script does the exact same as 'NerveSegmenetation3D.ipynb'
"""
import os
import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
from scipy.ndimage.interpolation import map_coordinates
from qimtools import visualization, inspection, io
from slgbuilder import GraphObject,MaxflowBuilder
def draw_segmentations(data, helper,layer=0):
"""Draw all segmentations for objects in the helper on top of the data."""
if data.ndim != 3:
raise ValueError('Data should be in three dimensions')
K = (len(helper.objects)//vol.shape[-1])//2
# Create figure.
plt.figure(figsize=(10, 10))
plt.imshow(data[...,layer], cmap='gray')
plt.xlim([0, data[...,layer].shape[1]-1])
plt.ylim([data[...,layer].shape[0]-1, 0])
# Draw segmentation lines.
for i, obj in enumerate(helper.objects):
if (i >= layer*K and i < (layer+1)*K) or (i >= (vol.shape[-1]+layer)*K and i< (vol.shape[-1]+layer+1)*K):
# 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()
def unfold_image(img, center, max_dists=None, r_min=1, r_max=20, angles=30, steps=15):
# 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
#%%
in_dir = 'data/'
Vol_path = os.path.join(in_dir, 'nerves3D.tiff')
# Load the data
vol = io.Volume( Vol_path )
# Convert the stack of 2D slices to a 3D volume
vol = vol.concatenate().astype(np.int32)
vol = np.transpose(vol,(1,2,0))
# visualization.show_vol( vol, axis=2 )
#%%
# 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 volume slice with centers.
plt.imshow(vol[...,0], 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)}')
#%%
nerve_samples = []
outer_nerves = []
inner_nerves = []
layers = []
# For each center, create an inner and outer never.
for i in range(vol.shape[-1]):
for center in centers_small:
samples, sample_points = unfold_image(vol[...,i], center,r_max=40,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(vol[...,i], 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))
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)
# Get segmentations.
segmentations = []
# for i in range(vol.shape[0]):
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]
# 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()
draw_segmentations(vol, helper,layer=0)
#%% Multi-object exclusion
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...')
# Add exclusion constraints between all pairs of outer nerves.
interval = len(outer_nerves)//vol.shape[-1]
for i in range(vol.shape[-1]):
for j in range(len(outer_nerves)//vol.shape[-1]):
for k in range(j + 1, len(outer_nerves)//vol.shape[-1]):
helper.add_layered_exclusion(outer_nerves[interval*i+j],
outer_nerves[interval*i+k], margin=3)
twice_flow = helper.solve()
print('Two times maximum flow/minimum energy:', twice_flow)
draw_segmentations(vol, helper,layer=9)
\ No newline at end of file
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