文章目录
待分割图片
直方图双峰阈值分割法
代码参考链接不需要进行tif->jpg，但保留处理代码，可用于处理高光谱数据。
代码：
import os
from osgeo import gdal
import numpy as np
import cv2
GRAY_SCALE = 256
# https://www.cnblogs.com/suoyike1001/p/15191940.html
def tif_jpg(rasterfile):
in_ds = gdal.Open(rasterfile) # 打开样本文件
xsize = in_ds.RasterXSize # 获取行列数
ysize = in_ds.RasterYSize
bands = in_ds.RasterCount
block_data = in_ds.ReadAsArray(0, 0, xsize, ysize).astype(np.float32)
B = block_data[0, :, :]
G = block_data[ 1,:, :]
R = block_data[2,:, :]
R1 = (R/np.max(R)*255)
G1 = (G / np.max(G) * 255)
B1 = (B / np.max(B) * 255)
data2 = cv2.merge([R1,G1,B1]).astype(np.int16)
return data2
def calcGrayHist(image):
'''
统计像素值
:param image:
:return:
'''
# 灰度图像的高，宽
rows, cols = image.shape
# 存储灰度直方图
grayHist = np.zeros([256], np.uint64)
for r in range(rows):
for c in range(cols):
grayHist[image[r][c]] += 1
return grayHist
#直方图全局阈值
def threshTwoPeaks(image):
# 计算灰度直方图
histogram = calcGrayHist(image)
# 找到灰度直方图的最大峰值对应的灰度值
maxLoc = np.where(histogram == np.max(histogram))
firstPeak = maxLoc[0][0]
# 寻找灰度直方图的第二个峰值对应的灰度值
measureDists = np.zeros([256], np.float32)
for k in range(256):
measureDists[k] = pow(k - firstPeak, 3) * histogram[k]#GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
maxLoc2 = np.where(measureDists == np.max(measureDists))
secondPeak = maxLoc2[0][0]
if firstPeak > secondPeak:
temp = histogram[int(secondPeak): int(firstPeak)]
minLoc = np.where(temp == np.min(temp))
thresh = secondPeak + minLoc[0][0]+ 1
else:
temp = histogram[int(firstPeak): int(secondPeak)]
minLoc = np.where(temp == np.min(temp))
thresh = firstPeak + minLoc[0][0]
img = image.copy()
img[img >= thresh] = 255
img[img < thresh] = 0
print("firstPeak:",firstPeak,",secondPeak:",secondPeak,",thresh:",thresh)
return img
def write_tiff(path,image_gray,out):
in_ds = gdal.Open(path) # 打开样本文件
xsize = in_ds.RasterXSize # 获取行列数
ysize = in_ds.RasterYSize
bands = in_ds.RasterCount
geotransform = in_ds.GetGeoTransform()
projection = in_ds.GetProjectionRef()
driver = gdal.GetDriverByName('GTiff')
outfile = out + "\\" + os.path.basename(path).split(".tif")[0] + "_mask.tif" # 对输出文件命名
out_dataset = driver.Create(outfile, xsize, ysize, 1, gdal.GDT_Float32) # 创建一个一波段的数据框架
out_band1 = out_dataset.GetRasterBand(1)
out_band1.WriteArray(image_gray)
out_dataset.SetGeoTransform(geotransform) # 写入仿射变换
out_dataset.SetProjection(projection)
return outfile
if __name__ == '__main__':
mask = r"190829vv.jpg"
out = r"ttt"
# img=tif_jpg(mask).astype(np.uint8)
img = cv2.imread(mask)
gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
seg_data = threshTwoPeaks(gray_img)
cv2.imwrite('histogram.jpg',seg_data)
# write_tiff(mask, seg_data , out)
分割结果：
阈值分割法
代码参考链接不需要进行tif->jpg，但保留处理代码，可用于处理高光谱数据。
代码：
import os
import cv2
import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal
GRAY_SCALE = 256
if __name__ == '__main__':
#设置阈值
thresh=40
mask = r"190811vv.jpg"
out = r"ttt"
# img=tif_jpg(mask).astype(np.uint8)
img = cv2.imread(mask)
gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
image_gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
#高斯滤波
img_blur = cv2.GaussianBlur(image_gray , (5, 5), 5)
# 阈值提取
img_blur[img_blur > thresh] = 255
img_blur[img_blur<= thresh] =1
# seg_data = threshTwoPeaks(gray_img)
cv2.imwrite('threshold.jpg',img_blur)
# write_tiff(mask, seg_data , out)
分割结果
模糊C均值
代码：
import random
import cv2
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['axes.unicode_minus']=False
plt.rcParams['font.sans-serif']=['SimHei']
import time
from PIL import Image
import torchvision.transforms as transforms
m=2 #隶属度函数因子
his_clust_num=[]
class FCM:
def __init__(self, data, clust_num,iter_num):
self.data = data
self.cnum = clust_num #聚类中心数
self.sample_num=data.shape[0] #样本集
self.dim = data.shape[-1] # 数据最后一维度数 特征数
Jlist=[] # 存储损失函数计算值的矩阵
U = self.Initial_U(self.sample_num, self.cnum) #初始化隶属度矩阵
for i in range(0, iter_num):
C = self.Cen_Iter(self.data, U, self.cnum)
U = self.U_Iter(U, C)
if (i==iter_num-1):
print("第%d次迭代" %(i+1))
print("聚类中心",C)
J = self.J_calcu(self.data, U, C) # 计算损失函数
Jlist = np.append(Jlist, J)
self.label = np.argmax(U, axis=0) # 所有样本的分类标签
self.Clast = C # 最后的类中心矩阵
self.Jlist = Jlist # 存储损失函数的矩阵
# 初始化隶属度矩阵U
def Initial_U(self, sample_num, cluster_n):
U = np.random.rand(sample_num, cluster_n) # sample_num为样本个数, cluster_n为分类数
row_sum = np.sum(U, axis=1) # 按行求和 row_sum: sample_num*1
row_sum = 1 / row_sum # 该矩阵每个数取倒数
U = np.multiply(U.T, row_sum) # 确保U的每列和为1 (cluster_n*sample_num).*(sample_num*1)
return U # cluster_n*sample_num
# 计算类中心
def Cen_Iter(self, data, U, cluster_n):
c_new = np.empty(shape=[0, self.dim]) # self.dim为样本矩阵的最后一维度
for i in range(0, cluster_n): # 图片像素值的dim为1
u_ij_m = U[i, :] ** m # (sample_num,)
sum_u = np.sum(u_ij_m)
ux = np.dot(u_ij_m, data) # (dim,)
ux = np.reshape(ux, (1, self.dim)) # (1,dim)
c_new = np.append(c_new, ux / sum_u, axis=0) # 按列的方向添加类中心到类中心矩阵
return c_new # cluster_num*dim
# 隶属度矩阵迭代
def U_Iter(self, U, c):
for i in range(0, self.cnum):
for j in range(0, self.sample_num):
sum = 0
for k in range(0, self.cnum):
temp = (np.linalg.norm(self.data[j, :] - c[i, :]) / np.linalg.norm(self.data[j, :] - c[k, :])) ** (2 / (m - 1))
sum = temp + sum
U[i, j] = 1 / sum
return U
# 计算损失函数值
def J_calcu(self, data, U, c):
temp1 = np.zeros(U.shape)
for i in range(0, U.shape[0]):
for j in range(0, U.shape[1]):
temp1[i, j] = (np.linalg.norm(data[j, :] - c[i, :])) ** 2 * U[i, j] ** m
J = np.sum(np.sum(temp1))
print("损失函数值:%.2f" %J)
return J
def FCM_pic_cut(img_path,gray=False,clustercenternum=3,iternum=5, outpath=None):
if gray:
img=cv2.imread(img_path,0) #灰度图
data=img.reshape(img.shape[0]*img.shape[1],1) #将图片拉成一列
else:
img=cv2.imread(img_path)
img=cv2.cvtColor(img,cv2.COLOR_BGR2RGB) #转化为RGB，不然plt时图片会怪怪的
data=img.reshape(-1,3) # 将三维降成二维
print("开始聚类")
start = time.time()
test=FCM(data,clustercenternum,iternum)
cluster=test.label # 聚类结果
center=test.Clast # 聚类中心
his_clust_num.append(test.Jlist[-1])
end = time.time()
print("聚类完成，开始生成图片")
#print(cluster)
new_img=center[cluster] # 根据聚类结果和聚类中心构建新图像
new_img=np.reshape(new_img,img.shape) #矩阵转成原来图片的形状
#print(new_img)
new_img=new_img.astype('uint8') # 要变成图像得数据得转换成uint8
plt.show()
new_img = Image.fromarray(np.uint8(new_img))
out_img = new_img.convert('L')
print('FPS:',1/(end - start))
plt.imsave(outpath,out_img) # 保存图片
if __name__ == '__main__':
# FCM_pic_cut("190811vv.jpg",gray=True,clustercenternum=3, outpath = '190811vvaoligti.jpg')
image_path = "190811vvaoligti.jpg"
# 读取图片
image = Image.open(image_path)
# 输出维度
print("RGB图像的维度：", np.array(image).shape)
# 显示原图
image.show()
# RGB转换我灰度图像
image_transforms = transforms.Compose([
transforms.Grayscale(1)
])
image = image_transforms(image)
# 输出灰度图像的维度
print("灰度图像维度： ", np.array(image).shape)
# 显示灰度图像
image.save("gray_fcm_190811vvaoligti.jpg.jpg")
输出：
开始聚类
损失函数值:168496691.92
损失函数值:168491447.17
损失函数值:168410643.84
损失函数值:167428097.59
第5次迭代
聚类中心 [[44.89842027]
[46.87565145]
[38.40223947]]
损失函数值:160078430.24
聚类完成，开始生成图片
FPS: 0.0018974307923545226
RGB图像的维度： (768, 1024, 3)
灰度图像维度： (768, 1024)
分割结果
超像素分割
代码链接
代码
# coding=utf-8
# 导入相应的python包
import argparse
from skimage import io
import matplotlib.pyplot as plt
from skimage.segmentation import slic
from skimage.util import img_as_float
from skimage.segmentation import mark_boundaries
# 设置并解析参数
ap = argparse.ArgumentParser()
ap.add_argument("-i", "--image", required = True, help = "Path to the image")
args = vars(ap.parse_args())
# 读取图片并将其转化为浮点型
image = img_as_float(io.imread(args["image"]))
# 循环设置不同的超像素组
for numSegments in (100, 200, 300):
# 应用slic算法并获取分割结果
segments = slic(image, n_segments = numSegments, sigma = 5)
# 绘制结果
fig = plt.figure("Superpixels -- %d segments" % (numSegments))
ax = fig.add_subplot(1, 1, 1)
ax.imshow(mark_boundaries(image, segments))
plt.axis("off")
plt.savefig("Superpixels%dsegments.jpg" % (numSegments))
# 显示结果
plt.show()
分割结果
K-means
代码链接
import numpy as np
def kmeans(image, num_clusters, seed=0, max_iter=10000):
np.random.seed(seed)
img = image.reshape(-1, 3)
index = 0
clus_value = np.array(np.random.rand(3, num_clusters)*255, dtype=float)
while True:
index += 1
flag = clus_value.copy()
cs = np.array([np.square(img-clus_value[0]), np.square(img-clus_value[1]), np.square(img-clus_value[2])])
cs = np.sum(cs, axis=2)
labels = np.argmin(cs, axis=0)
for j in range(num_clusters):
clus_value[j] = img[labels==j].mean(axis=0)
if np.sum(np.abs(clus_value-flag)) < 1e-8 or index==max_iter:
break
segmented_image = np.array(clus_value[labels].reshape(image.shape), dtype=np.uint8)
# print(segmented_image.shape,np.argmin(segmented_image))
segmented_image[segmented_image>np.mean(segmented_image)] = 255
return segmented_image
import cv2
img = cv2.imread("190811vv.jpg")
segment_image = kmeans(img, num_clusters=3)
cv2.imwrite("filename11.png", segment_image)
分割结果
主动轮廓模型
snake代码
# -*- coding: utf-8 -*-
"""
====================
Morphological Snakes
====================
代码来源
https://blog.csdn.net/lironglu_/article/details/123917371
"""
__author__ = "P. Márquez Neila <p.mneila@upm.es>"
from itertools import cycle
import numpy as np
from scipy import ndimage as ndi
__all__ = [
'morphological_chan_vese',
'morphological_geodesic_active_contour',
'inverse_gaussian_gradient',
'circle_level_set',
'checkerboard_level_set'
]
__version__ = (2, 1, 1)
__version_str__ = ".".join(map(str, __version__))
class _fcycle(object):
def __init__(self, iterable):
"""Call functions from the iterable each time it is called."""
self.funcs = cycle(iterable)
def __call__(self, *args, **kwargs):
f = next(self.funcs)
return f(*args, **kwargs)
# SI and IS operators for 2D and 3D.
_P2 = [np.eye(3),
np.array([[0, 1, 0]] * 3, dtype=object),
np.flipud(np.eye(3)),
np.rot90([[0, 1, 0]] * 3)]
_P3 = [np.zeros((3, 3, 3), dtype=object) for i in range(9)]
_P3[0][:, :, 1] = 1
_P3[1][:, 1, :] = 1
_P3[2][1, :, :] = 1
_P3[3][:, [0, 1, 2], [0, 1, 2]] = 1
_P3[4][:, [0, 1, 2], [2, 1, 0]] = 1
_P3[5][[0, 1, 2], :, [0, 1, 2]] = 1
_P3[6][[0, 1, 2], :, [2, 1, 0]] = 1
_P3[7][[0, 1, 2], [0, 1, 2], :] = 1
_P3[8][[0, 1, 2], [2, 1, 0], :] = 1
def sup_inf(u):
"""SI operator."""
if np.ndim(u) == 2:
P = _P2
elif np.ndim(u) == 3:
P = _P3
else:
raise ValueError("u has an invalid number of dimensions "
"(should be 2 or 3)")
erosions = []
for P_i in P:
erosions.append(ndi.binary_erosion(u, P_i))
return np.array(erosions, dtype=np.int8).max(0)
def inf_sup(u):
"""IS operator."""
if np.ndim(u) == 2:
P = _P2
elif np.ndim(u) == 3:
P = _P3
else:
raise ValueError("u has an invalid number of dimensions "
"(should be 2 or 3)")
dilations = []
for P_i in P:
dilations.append(ndi.binary_dilation(u, P_i))
return np.array(dilations, dtype=np.int8).min(0)
_curvop = _fcycle([lambda u: sup_inf(inf_sup(u)), # SIoIS
lambda u: inf_sup(sup_inf(u))]) # ISoSI
def _check_input(image, init_level_set):
"""Check that shapes of `image` and `init_level_set` match."""
if image.ndim not in [2, 3]:
raise ValueError("`image` must be a 2 or 3-dimensional array.")
if len(image.shape) != len(init_level_set.shape):
raise ValueError("The dimensions of the initial level set do not "
"match the dimensions of the image.")
def _init_level_set(init_level_set, image_shape):
if isinstance(init_level_set, str):
if init_level_set == 'checkerboard':
res = checkerboard_level_set(image_shape)
elif init_level_set == 'circle':
res = circle_level_set(image_shape)
elif init_level_set == 'ellipsoid':
res = ellipsoid_level_set(image_shape)
else:
raise ValueError("`init_level_set` not in "
"['checkerboard', 'circle', 'ellipsoid']")
else:
res = init_level_set
return res
def circle_level_set(image_shape, center=None, radius=None):
if center is None:
center = tuple(i // 2 for i in image_shape)
if radius is None:
radius = min(image_shape) * 3.0 / 8.0
grid = np.mgrid[[slice(i) for i in image_shape]]
grid = (grid.T - center).T
phi = radius - np.sqrt(np.sum((grid)**2, 0))
res = np.int8(phi > 0)
return res
def ellipsoid_level_set(image_shape, center=None, semi_axis=None):
if center is None:
center = tuple(i // 2 for i in image_shape)
if semi_axis is None:
semi_axis = tuple(i / 2 for i in image_shape)
if len(center) != len(image_shape):
raise ValueError("`center` and `image_shape` must have the same length.")
if len(semi_axis) != len(image_shape):
raise ValueError("`semi_axis` and `image_shape` must have the same length.")
if len(image_shape) == 2:
xc, yc = center
rx, ry = semi_axis
phi = 1 - np.fromfunction(
lambda x, y: ((x - xc) / rx) ** 2 +
((y - yc) / ry) ** 2,
image_shape, dtype=float)
elif len(image_shape) == 3:
xc, yc, zc = center
rx, ry, rz = semi_axis
phi = 1 - np.fromfunction(
lambda x, y, z: ((x - xc) / rx) ** 2 +
((y - yc) / ry) ** 2 +
((z - zc) / rz) ** 2,
image_shape, dtype=float)
else:
raise ValueError("`image_shape` must be a 2- or 3-tuple.")
res = np.int8(phi > 0)
return res
def checkerboard_level_set(image_shape, square_size=5):
grid = np.ogrid[[slice(i) for i in image_shape]]
grid = [(grid_i // square_size) & 1 for grid_i in grid]
checkerboard = np.bitwise_xor.reduce(grid, axis=0, )
res = np.int8(checkerboard)
return res
def inverse_gaussian_gradient(image, alpha=100.0, sigma=5.0):
gradnorm = ndi.gaussian_gradient_magnitude(image, sigma, mode='nearest')
return 1.0 / np.sqrt(1.0 + alpha * gradnorm)
def morphological_chan_vese(image, iterations, init_level_set='checkerboard',
smoothing=1, lambda1=1, lambda2=1,
iter_callback=lambda x: None):
init_level_set = _init_level_set(init_level_set, image.shape)
_check_input(image, init_level_set)
u = np.int8(init_level_set > 0)
iter_callback(u)
for _ in range(iterations):
# inside = u > 0
# outside = u <= 0
c0 = (image * (1 - u)).sum() / float((1 - u).sum() + 1e-8)
c1 = (image * u).sum() / float(u.sum() + 1e-8)
# Image attachment
du = np.gradient(u)
abs_du = np.abs(du).sum(0)
aux = abs_du * (lambda1 * (image - c1)**2 - lambda2 * (image - c0)**2)
u[aux < 0] = 1
u[aux > 0] = 0
# Smoothing
for _ in range(smoothing):
u = _curvop(u)
iter_callback(u)
return u
def morphological_geodesic_active_contour(gimage , iterations,
init_level_set='circle', smoothing=1,
threshold='auto', balloon=0,
iter_callback=lambda x: None):
image = gimage
init_level_set = _init_level_set(init_level_set, image.shape)
_check_input(image, init_level_set)
if threshold == 'auto':
threshold = np.percentile(image, 40)
structure = np.ones((3,) * len(image.shape), dtype=np.int8)
dimage = np.gradient(image)
# threshold_mask = image > threshold
if balloon != 0:
threshold_mask_balloon = image > threshold / np.abs(balloon)
u = np.int8(init_level_set > 0)
iter_callback(u)
for _ in range(iterations):
# Balloon
if balloon > 0:
aux = ndi.binary_dilation(u, structure)
elif balloon < 0:
aux = ndi.binary_erosion(u, structure)
if balloon != 0:
u[threshold_mask_balloon] = aux[threshold_mask_balloon]
# Image attachment
aux = np.zeros_like(image)
du = np.gradient(u)
for el1, el2 in zip(dimage, du):
aux += el1 * el2
u[aux > 0] = 1
u[aux < 0] = 0
# Smoothing
for _ in range(smoothing):
u = _curvop(u)
iter_callback(u)
return u
def No_morphological_geodesic_active_contour(gimage , iterations,
init_level_set='circle', smoothing=1,
threshold='auto', balloon=0):
image = gimage
init_level_set = _init_level_set(init_level_set, image.shape)
_check_input(image, init_level_set)
if threshold == 'auto':
threshold = np.percentile(image, 40)
structure = np.ones((3,) * len(image.shape), dtype=np.int8)
dimage = np.gradient(image)
# threshold_mask = image > threshold
if balloon != 0:
threshold_mask_balloon = image > threshold / np.abs(balloon)
u = np.int8(init_level_set > 0)
for _ in range(iterations):
# Balloon
if balloon > 0:
aux = ndi.binary_dilation(u, structure)
elif balloon < 0:
aux = ndi.binary_erosion(u, structure)
if balloon != 0:
u[threshold_mask_balloon] = aux[threshold_mask_balloon]
# Image attachment
aux = np.zeros_like(image)
du = np.gradient(u)
for el1, el2 in zip(dimage, du):
aux += el1 * el2
u[aux > 0] = 1
u[aux < 0] = 0
# Smoothing
for _ in range(smoothing):
u = _curvop(u)
return u
利用canny提取图片边缘
import cv2
import numpy as np
import logging
from imageio import imread
import matplotlib
from matplotlib import pyplot as plt
import morphsnakes as ms # 调用morphsnakes.py
image = cv2.imread('190811vv.jpg')# 输入图片名称,路径
img = image
threshold = np.percentile(image, 40)# 设定一个基于40% 的 阈值
canny_img = cv2.Canny(img,threshold,255)
cv2.imshow('canny',canny_img)
cv2.imwrite('canny.jpg',canny_img)
cv2.waitKey(0)
对图片进行闭运算, 获得更明确的边缘信息
# 闭运算
def close_img(img,k):
if len(img.shape)==3:
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
ret, binary = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY_INV | cv2.THRESH_OTSU)
else:
binary = img
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (k, k))
binary = cv2.morphologyEx(binary, cv2.MORPH_CLOSE, kernel)
return binary
################################################################
kernel_size = 9
img = canny_img
closed_img = close_img(img,kernel_size)
cv2.imwrite('closed_img.jpg', closed_img)
cv2.imshow('closed_img', closed_img)
cv2.waitKey(0)
采用snake算法收缩
def snake(img, iterations = 2000):
gimg = ms.inverse_gaussian_gradient(img, alpha=500, sigma=4.0) # 反高斯梯度
# 初始轮廓设置
init_ls = np.zeros(img.shape, dtype=np.int8)
init_ls[10:-10, 10:-10] = 1
# MorphGAC.
uuu = ms.No_morphological_geodesic_active_contour(gimg, iterations, init_ls,
smoothing=5, threshold=0.1,
balloon=-1)# 主动轮廓模型输出是0-1矩阵, 需要进行灰度线性拉伸
uuu = grey_scale(uuu)
return uuu
# 灰度拉伸
def grey_scale(img):
rows, cols = img.shape
flat_gray = img.reshape((cols * rows,)).tolist()
min1 = min(flat_gray)
max1 = max(flat_gray)
# print('min = %d,max = %d' % (min1, max1))
output = np.uint8(255 / (max1 - min1) * (img - min1) + 0.5)
return output
seged_img = snake(closed_img, iterations = 1000)
cv2.imshow('seged',seged_img)
cv2.waitKey(0)
效果不好, 需要调整参数, 但实际上在闭运算处就已经完成部分分割任务了, 这个实际上可以在闭运算kernel比较小, 有空心区域的时候用主动轮廓模型分割效果会好一些
667
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
