正确关闭使用 skimage.measure.find_contours() 生成的多边形
Correctly closing of polygons generated using skimage.measure.find_contours()
我目前正在使用 skimage.measure.find_contours()
在表面上查找轮廓。现在我已经找到了轮廓,我需要能够找到它们所包含的区域。
当所有顶点都在数据集中时,这很好,因为具有完全封闭的多边形。
但是,如果轮廓在边缘或角处突破表面的边缘,我如何确保多边形完全封闭?发生这种情况时,我想使用曲面的边缘作为附加顶点来关闭多边形。例如在下图中,显示了等高线,您可以看到等高线在图像的边缘结束,我该如何关闭它们?同样在棕色轮廓的例子中,它只是一条线,我不认为我想要返回一个区域,我如何挑出这种情况?
我知道我可以通过检查多边形的最后一个顶点是否与第一个顶点相同来检查封闭的 contours/polygons。
我有计算多边形内面积的代码,取自here
def find_area(array):
a = 0
ox,oy = array[0]
for x,y in array[1:]:
a += (x*oy-y*ox)
ox,oy = x,y
return -a/2
我只需要关闭多边形的帮助。并检查可能发生的不同情况。
谢谢
更新:
应用@soupault 建议的解决方案后,我得到了这个代码:
import numpy as np
import matplotlib.pyplot as plt
from skimage import measure
# Construct some test data
x, y = np.ogrid[-np.pi:np.pi:100j, -np.pi:np.pi:100j]
r = np.sin(np.exp((np.sin(x)**3 + np.cos(y)**2)))
# Coordinates of point of interest
pt = [(49,75)]
# Apply thresholding to the surface
threshold = 0.8
blobs = r > threshold
# Make a labelled image based on the thresholding regions
blobs_labels = measure.label(blobs, background = 0)
# Show the thresholded regions
plt.figure()
plt.imshow(blobs_labels, cmap='spectral')
# Apply regionprops to charactersie each of the regions
props = measure.regionprops(blobs_labels, intensity_image = r)
# Loop through each region in regionprops, identify if the point of interest is
# in that region. If so, plot the region and print it's area.
plt.figure()
plt.imshow(r, cmap='Greys')
plt.plot(pt[0][0], pt[0][1],'rx')
for prop in props:
coords = prop.coords
if np.sum(np.all(coords[:,[1,0]] == pt[0], axis=1)):
plt.plot(coords[:,1],coords[:,0],'r.')
print(prop.area)
此解决方案假定每个像素的大小为 1x1。在我的真实数据解决方案中,情况并非如此,因此我还应用了以下函数对数据应用线性插值。我相信你也可以应用类似的功能,使每个像素的面积更小,从而提高数据的分辨率。
import numpy as np
from scipy import interpolate
def interpolate_patch(x,y,patch):
x_interp = np.arange(np.ceil(x[0]), x[-1], 1)
y_interp = np.arange(np.ceil(y[0]), y[-1], 1)
f = interpolate.interp2d(x, y, patch, kind='linear')
patch_interp = f(x_interp, y_interp)
return x_interp, y_interp, patch_interp
如果您需要测量不同区域的属性,自然会从找到区域(而不是轮廓)开始。
算法如下,在本例中:
准备一张带标签的图片:
1.a或者用不同颜色填充不同轮廓线之间的区域;
1.b 或者应用一些图像阈值函数,然后 运行 skimage.measure.label
(http://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.label);
执行 regionprops
使用非常标记的图像作为输入 (http://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.regionprops);
遍历 regionprops
中的区域并计算所需的参数(面积、周长等)。
通过 regionprops
识别图像中的区域后,您可以为每个区域调用 .coords
以获得封闭的轮廓。
如果有人需要通过图像边缘关闭开放轮廓(并制作多边形),这里是:
import shapely.geometry as sgeo
import shapely.ops as sops
def close_contour_with_image_edge(contour, image_shape):
"""
this function uses shapely because its easiest way to do that
:param contour: contour generated by skimage.measure.find_contours()
:param image_shape: tuple (row, cols), standard return of numpy shape()
:return:
"""
# make contour linestring
contour_line = sgeo.LineString(contour)
# make image box linestring
box_rows, box_cols = image_shape[0], image_shape[1]
img_box = sgeo.LineString(coordinates=(
(0, 0),
(0, box_cols-1),
(box_rows-1, box_cols-1),
(box_rows-1, 0),
(0, 0)
))
# intersect box with non-closed contour and get shortest line which touch both of contour ends
edge_points = img_box.intersection(contour_line)
edge_parts = sops.split(img_box, edge_points)
edge_parts = list(part for part in edge_parts.geoms if part.touches(edge_points.geoms[0]) and part.touches(edge_points.geoms[1]))
edge_parts.sort(reverse=False, key=lambda x: x.length)
contour_edge = edge_parts[0]
# weld it
contour_line = contour_line.union(contour_edge)
contour_line = sops.linemerge(contour_line)
contour_polygon = sgeo.Polygon(contour_line.coords)
return contour_polygon
我目前正在使用 skimage.measure.find_contours()
在表面上查找轮廓。现在我已经找到了轮廓,我需要能够找到它们所包含的区域。
当所有顶点都在数据集中时,这很好,因为具有完全封闭的多边形。 但是,如果轮廓在边缘或角处突破表面的边缘,我如何确保多边形完全封闭?发生这种情况时,我想使用曲面的边缘作为附加顶点来关闭多边形。例如在下图中,显示了等高线,您可以看到等高线在图像的边缘结束,我该如何关闭它们?同样在棕色轮廓的例子中,它只是一条线,我不认为我想要返回一个区域,我如何挑出这种情况?
我知道我可以通过检查多边形的最后一个顶点是否与第一个顶点相同来检查封闭的 contours/polygons。
我有计算多边形内面积的代码,取自here
def find_area(array):
a = 0
ox,oy = array[0]
for x,y in array[1:]:
a += (x*oy-y*ox)
ox,oy = x,y
return -a/2
我只需要关闭多边形的帮助。并检查可能发生的不同情况。
谢谢
更新:
应用@soupault 建议的解决方案后,我得到了这个代码:
import numpy as np
import matplotlib.pyplot as plt
from skimage import measure
# Construct some test data
x, y = np.ogrid[-np.pi:np.pi:100j, -np.pi:np.pi:100j]
r = np.sin(np.exp((np.sin(x)**3 + np.cos(y)**2)))
# Coordinates of point of interest
pt = [(49,75)]
# Apply thresholding to the surface
threshold = 0.8
blobs = r > threshold
# Make a labelled image based on the thresholding regions
blobs_labels = measure.label(blobs, background = 0)
# Show the thresholded regions
plt.figure()
plt.imshow(blobs_labels, cmap='spectral')
# Apply regionprops to charactersie each of the regions
props = measure.regionprops(blobs_labels, intensity_image = r)
# Loop through each region in regionprops, identify if the point of interest is
# in that region. If so, plot the region and print it's area.
plt.figure()
plt.imshow(r, cmap='Greys')
plt.plot(pt[0][0], pt[0][1],'rx')
for prop in props:
coords = prop.coords
if np.sum(np.all(coords[:,[1,0]] == pt[0], axis=1)):
plt.plot(coords[:,1],coords[:,0],'r.')
print(prop.area)
此解决方案假定每个像素的大小为 1x1。在我的真实数据解决方案中,情况并非如此,因此我还应用了以下函数对数据应用线性插值。我相信你也可以应用类似的功能,使每个像素的面积更小,从而提高数据的分辨率。
import numpy as np
from scipy import interpolate
def interpolate_patch(x,y,patch):
x_interp = np.arange(np.ceil(x[0]), x[-1], 1)
y_interp = np.arange(np.ceil(y[0]), y[-1], 1)
f = interpolate.interp2d(x, y, patch, kind='linear')
patch_interp = f(x_interp, y_interp)
return x_interp, y_interp, patch_interp
如果您需要测量不同区域的属性,自然会从找到区域(而不是轮廓)开始。
算法如下,在本例中:
准备一张带标签的图片:
1.a或者用不同颜色填充不同轮廓线之间的区域;
1.b 或者应用一些图像阈值函数,然后 运行
skimage.measure.label
(http://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.label);执行
regionprops
使用非常标记的图像作为输入 (http://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.regionprops);遍历
regionprops
中的区域并计算所需的参数(面积、周长等)。
通过 regionprops
识别图像中的区域后,您可以为每个区域调用 .coords
以获得封闭的轮廓。
如果有人需要通过图像边缘关闭开放轮廓(并制作多边形),这里是:
import shapely.geometry as sgeo
import shapely.ops as sops
def close_contour_with_image_edge(contour, image_shape):
"""
this function uses shapely because its easiest way to do that
:param contour: contour generated by skimage.measure.find_contours()
:param image_shape: tuple (row, cols), standard return of numpy shape()
:return:
"""
# make contour linestring
contour_line = sgeo.LineString(contour)
# make image box linestring
box_rows, box_cols = image_shape[0], image_shape[1]
img_box = sgeo.LineString(coordinates=(
(0, 0),
(0, box_cols-1),
(box_rows-1, box_cols-1),
(box_rows-1, 0),
(0, 0)
))
# intersect box with non-closed contour and get shortest line which touch both of contour ends
edge_points = img_box.intersection(contour_line)
edge_parts = sops.split(img_box, edge_points)
edge_parts = list(part for part in edge_parts.geoms if part.touches(edge_points.geoms[0]) and part.touches(edge_points.geoms[1]))
edge_parts.sort(reverse=False, key=lambda x: x.length)
contour_edge = edge_parts[0]
# weld it
contour_line = contour_line.union(contour_edge)
contour_line = sops.linemerge(contour_line)
contour_polygon = sgeo.Polygon(contour_line.coords)
return contour_polygon