In [34]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import imageio
%matplotlib inline
from PIL import Image
from scipy import ndimage as ndi
from skimage.morphology import watershed, binary_closing, binary_opening, binary_erosion, binary_dilation, erosion
from skimage.feature import peak_local_max
from skimage.filters import threshold_otsu, sobel, laplace, gaussian, median
In [73]:
def list_plot(im_list, title_list, figsize=(20, 10), cmap=None, ax_on=False):
# Plot a list of image with titles
fig, axis = plt.subplots(1, len(im_list), figsize=figsize)
for i in range(len(im_list)):
if isinstance(cmap, list):
c = cmap[i]
else:
c = cmap
axis[i].imshow(im_list[i], cmap=c)
axis[i].set_title(title_list[i])
if not ax_on:
axis[i].set_axis_off()
def plot_histogram(f, nbins, title='', color='b'):
h, bin_edges = np.histogram(f, nbins,(0,180))
# Plot
w=255./nbins
bin_centers = bin_edges[1:]-(w/2)
plt.figure()
plt.bar(bin_centers, h, width=w, color=color)
plt.title(title)
# normalize images to specified range
def normalize(f, range=(0,255)):
# normalize to 0,1 first
f = (f - np.min(f)) / (np.max(f) - np.min(f))
# set to specified range
f = (f + np.min(range)) * (np.max(range) - np.min(range))
return f
def shuffle_labels(labels):
u = np.unique(labels)
idx = labels == 0
np.random.shuffle(u)
out = u[labels] + 1
out[idx] = 0
return out
def make_gif(images, path, add_color=True, cmap='gnuplot'):
cm = matplotlib.cm.get_cmap(cmap)
if add_color:
images = [cm(normalize(im).astype(int)) for im in images]
imageio.mimsave(path, images)
A implementação do algorítimo watershed sem marcadores ocorre iterativamente inundando os pixels de menor valor da imagem até que todos os pixels sejam inundados. Caso ocorra a mistura de dois catchment basins (rótulos diferentes) uma fronteira/barreira é criada entre eles usando operações morfológicas (dilatação).
A cada iteração utiliza-se a função label do pacote scipy para identificar os catchment basin (regiões conexas na imagem). Em seguida, calcula-se a intersecção de cada elemento $q$ encontrado na iteração $N$ com a imagem inundada da iteração anterior $N-1$. A intersecção, por sua vez, também pode ser rotulada em elementos contínuos e total de elementos $c$ definem o próximo passo do algorítimo.
catchement basin. Ou seja $q$ que recebe um novo rotulo.
In [3]:
def custom_watershed(image, verbose=False, return_intermediate=False):
im = normalize(image.copy()).astype(int)
unique_values = np.unique(im)
kernel = np.ones((5,5))
structure = np.ones((3,3))
out = []
for i in range(len(unique_values)):
if not i:
continue
labels = ndi.label(im < unique_values[i], structure)[0]
if return_intermediate:
out.append(labels)
prev_step = im < unique_values[i-1]
for l in np.unique(labels):
if not l:
continue
q = labels == l
overlap = ndi.label(q & prev_step, structure)[0]
u_overlap = np.unique(overlap)
count = len(u_overlap) - 1
if count == 0:
if verbose:
print('New object found')
elif count == 1:
if verbose:
print('Increasing region')
elif count > 1:
if verbose:
print('Objects merged, constructing dam')
dam = np.zeros_like(im)
for d in u_overlap:
if not d:
continue
qi = overlap == d
dam += binary_dilation(qi, kernel)
im[dam > 1] = np.max(im) + 1
if return_intermediate:
return ndi.label(im < 255)[0], out
else:
return ndi.label(im < 255)[0]
In [111]:
x, y = np.indices((256, 256))
x1, y1 = 89, 89
x2, y2 = 140, 166
x3, y3 = 54, 160
r1, r2, r3 = 48, 64, 32
mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2
mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2
mask_circle3 = (x - x3)**2 + (y - y3)**2 < r3**2
image = np.logical_or(mask_circle1, mask_circle2)
image = np.logical_or(image, mask_circle3)
plt.imshow(image, cmap='gray')
plt.title('Image sintética criada')
Out[111]:
Para transformar a imagem binária em uma imagem em tons de cinzas análoga a um relevo, utiliza-se a transformada de distância.
Nota-se que tal transformada forma morros e picos ou invés de vales (basins), desta forma utiliza-se a imagem transformada multiplicada por -1.
In [112]:
distance = ndi.distance_transform_edt(image)
list_plot([distance, -distance], 'Distância_Distância invertida'.split('_'), cmap='gray', figsize=(10,5))
In [113]:
seg_im = custom_watershed(-distance)
plt.imshow(seg_im, cmap='gnuplot')
Out[113]:
In [69]:
gray_im = np.asarray(Image.open('metal.jpg').convert('L'))
plt.imshow(gray_im, cmap='gray')
plt.title('Imagem RGB')
Out[69]:
In [70]:
T = 150
bin_im = gray_im > T
plt.imshow(bin_im, cmap='gray')
plt.title('Imagem binária')
Out[70]:
In [71]:
distance = ndi.distance_transform_edt(bin_im)
list_plot([distance, -distance],
'Resultado da transformada_Resultado da transformada multiplicado por -1'.split('_'),
cmap='gray', figsize=(15, 7))
In [77]:
seg_im = custom_watershed(-distance)
plt.figure(figsize=(10, 10))
plt.imshow(shuffle_labels(seg_im), cmap='jet')
print('{} objetos encontrados.' .format(len(np.unique(seg_im))))
Primeiro abre-se a foto original e converte-se para HSV (será aplicado uma limiarização no canal de saturação)
In [78]:
image = Image.open('moedas.jpg')
hsv_im = np.asarray(image.convert('HSV'))
gray_im = np.asarray(image.convert('L'))
plt.imshow(image)
plt.title('Imagem RGB')
Out[78]:
Visualiza-se os histogramas dos três canais da imagem.
In [79]:
plot_histogram(hsv_im[:, :, 0], 40, 'Histograma da matiz', color='b')
plot_histogram(hsv_im[:, :, 1], 40, 'Histograma da saturação', color='r')
plot_histogram(hsv_im[:, :, 2], 40, 'Histograma do valor', color='g')
Nota-se claramente uma distribuição quase bimodal no histograma da saturação. Pode-se aproveitar tal fato para realizar uma limiarização na imagem a partir deste canal.
In [80]:
threshold = 40
bin_im = hsv_im[:, :, 1] > threshold
plt.imshow(bin_im, cmap='gray')
plt.title('Imagem binária depois da limiarização')
Out[80]:
In [83]:
distance = ndi.distance_transform_edt(bin_im)
seg_im = custom_watershed(-distance)
plt.imshow(shuffle_labels(seg_im), cmap='jet')
Out[83]:
Nota-se que ocorreu um processo de oversegmentation. Isto ocorre devido a minimos locais existentes em uma imagem real.
In [18]:
x, y = np.indices((80, 80))
x1, y1 = 28, 28
x2, y2 = 44, 52
x3, y3 = 17, 50
r1, r2, r3 = 16, 20, 10
mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2
mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2
mask_circle3 = (x - x3)**2 + (y - y3)**2 < r3**2
image = np.logical_or(mask_circle1, mask_circle2)
image = np.logical_or(image, mask_circle3)
plt.imshow(image, cmap='gray')
plt.title('Image sintética criada')
Out[18]:
In [19]:
distance = ndi.distance_transform_edt(image)
local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((5, 5)),
labels=image)
im_list = [image, distance, -distance, local_maxi]
title_list = 'Imagem original_Distâncias_Distâncias * -1_Máximo local (distâncias)'.split('_')
list_plot(im_list, title_list, cmap='gray')
In [20]:
markers = ndi.label(local_maxi)[0]
seg_im = watershed(-distance, markers, mask=image)
im_list = [image, distance, markers, seg_im]
title_list = 'Imagem original_Distâncias_Marcadores_Imagem segmentada'.split('_')
cmap_list = ['gray', 'gray', 'gnuplot', 'gnuplot']
list_plot(im_list, title_list, cmap=cmap_list)
In [21]:
gray_im = np.asarray(Image.open('metal.jpg').convert('L'))
plt.imshow(gray_im, cmap='gray')
plt.title('Imagem RGB')
Out[21]:
In [22]:
T = 150
bin_im = gray_im > T
plt.imshow(bin_im, cmap='gray')
plt.title('Imagem binária')
Out[22]:
Agora usa-se a transformada de distâncias para criar a imagem relevo a partir da imagem binária. Apartir do resultado, o algoritimo de busca de picos locais é usado para criar os marcadores. Finalmente os objetos são segmentados usando o algorítimo de watershed no inverso das distâncias com auxilio dos marcadores.
In [23]:
distance = ndi.distance_transform_edt(bin_im)
kernel = 10
local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((kernel, kernel)))
markers = ndi.label(local_maxi)[0]
seg_im = watershed(-distance, markers, mask=bin_im)
im_list = [gray_im, distance, shuffle_labels(seg_im)]
title_list = 'Imagem original_Distâncias_Imagem segmentada'.split('_')
cmap_list = ['gray', 'gray', 'jet', 'jet']
list_plot(im_list, title_list, figsize=(20, 20), cmap=cmap_list)
print('{} objetos localizados.'.format(len(np.unique(seg_im))) )
Primeiro abre-se a foto original e converte-se para HSV (será aplicado uma limiarização no canal de saturação)
In [24]:
image = Image.open('moedas.jpg')
hsv_im = np.asarray(image.convert('HSV'))
gray_im = np.asarray(image.convert('L'))
plt.imshow(image)
plt.title('Imagem RGB')
Out[24]:
Visualiza-se os histogramas dos três canais da imagem.
In [25]:
plot_histogram(hsv_im[:, :, 0], 40, 'Histograma da matiz', color='b')
plot_histogram(hsv_im[:, :, 1], 40, 'Histograma da saturação', color='r')
plot_histogram(hsv_im[:, :, 2], 40, 'Histograma do valor', color='g')
Nota-se claramente uma distribuição quase bimodal no histograma da saturação. Pode-se aproveitar tal fato para realizar uma limiarização na imagem a partir deste canal.
In [26]:
threshold = 40
bin_im = hsv_im[:, :, 1] > threshold
plt.imshow(bin_im, cmap='gray')
plt.title('Imagem binária depois da limiarização')
Out[26]:
Agora usa-se a transformada de distâncias para criar a imagem relevo a partir da imagem binária. Apartir do resultado, o algoritimo de busca de picos locais é usado para criar os marcadores. Finalmente os objetos são segmentados usando o algorítimo de watershed no inverso das distâncias com auxilio dos marcadores.
In [27]:
distance = ndi.distance_transform_edt(bin_im)
kernel = 15
local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((kernel, kernel)))
markers = ndi.label(local_maxi)[0]
seg_im = watershed(-distance, markers, mask=bin_im)
im_list = [gray_im, distance, seg_im]
title_list = 'Imagem original_Distâncias_Imagem segmentada'.split('_')
cmap_list = ['gray', 'gray', 'jet', 'jet']
list_plot(im_list, title_list, figsize=(20, 20), cmap=cmap_list)
Contando os valores únicos da imagem segmentada é possível contar o número de objetos encotrados
In [28]:
print('{} objetos identificados de um total de 32.'.format((np.unique(seg_im) > 0).sum()))
plt.imshow(gray_im, cmap='gray')
axs = plt.gca()
seg_im = shuffle_labels(seg_im).astype(float)
seg_im[seg_im == 0] = np.nan
p = axs.imshow(seg_im, cmap='jet', alpha=0.5)
plt.colorbar(p)
Out[28]: