In [1]:
%pylab inline
from sympy import init_printing
init_printing(use_latex='mathjax')
import os
os.chdir(os.path.join(os.curdir,'data'))
In [2]:
from PIL import Image
pil_im = Image.open('empire.jpg')
pil_im
Out[2]:
In [3]:
import cv2
cv_im = cv2.imread('empire.jpg')
cv_im = cv2.cvtColor(cv_im,cv2.COLOR_BGR2RGB) # 转为RGB用于matplot显示
imshow(cv_im)
Out[3]:
In [4]:
Image.open('empire.jpg').convert('L')
Out[4]:
In [5]:
thumbnail = pil_im.copy()
thumbnail.thumbnail((128,128)) # 缩略图
thumbnail
Out[5]:
In [6]:
box = (100,100,400,400) # (左,上,右,下)
region = pil_im.crop(box) # 剪裁指定区域
region = region.transpose(Image.ROTATE_180) # 转180度
tmp = pil_im.copy()
tmp.paste(region,box) # 黏贴回去
tmp
Out[6]:
In [7]:
pil_im.resize((128,128)) # 生成指定大小的新图像
Out[7]:
In [8]:
pil_im.rotate(45).resize((128,128)) # 逆时针转45度
Out[8]:
In [9]:
%pylab inline
In [10]:
from PIL import Image
from pylab import *
im = array(Image.open('empire.jpg'))
imshow(im)
x = [100, 100, 400, 400]
y = [200, 500, 200, 500]
plot(x, y, 'r*')
plot(x[:2], y[:2])
#axis('off')
title('Plotting: "empire.jpg"')
show()
In [11]:
from PIL import Image
from pylab import *
im = array(Image.open('empire.jpg'))
imshow(im)
x = [100, 100, 400, 400]
y = [200, 500, 200, 500]
plot(x, y, 'r*')
plot(x[:2], y[:2])
axis('off')
title('Plotting: "empire.jpg"')
show()
绘图格式命令:
w : white 白
- : 实线
: : 点线
. : 点
In [12]:
im = pil_im.convert('L')
figure()
gray() # 指定为黑白色
contour(im, orign='image') # 轮廓图像
axis('equal') # 不变形
axis('off') # 不显示
Out[12]:
In [13]:
figure()
np_im = np.array(im)
hist(np_im.flatten(),128) # 直方图
xlim([0,256])
show()
In [14]:
# cv2.calcHist(images, channels, mask, histSize, ranges[, hist[, accumulate ]]) -> hist
import cv2
im = cv2.cvtColor(array(pil_im),cv2.COLOR_BGR2GRAY) # 转为灰度图像
hist_item = cv2.calcHist([im],[0],None,[128],[0,256]) # 计算直方图
xlim([0,256])
bar(array([x for x in range(128)])*2,hist_item)
Out[14]:
In [15]:
hist_item = cv2.calcHist([im],[0],None,[256],[0,256]) # 计算直方图
cv2.normalize(hist_item,hist_item,0,255,cv2.NORM_MINMAX) # 归一化到 [0,255]的范围
hist_item=np.int32(np.around(hist_item)) # 像素取整
ylim([0,256])
xlim([0,256])
bar([x for x in range(256)],hist_item)
Out[15]:
In [16]:
## 交互式标注
# im = pil_im.copy()
# x = ginput(3) # 点击3次选择3个点
# x
In [17]:
im = np.array(pil_im.copy())
i,j,k = 0,0,0
im[i,j,k] # 坐标 i行,j列,颜色通道 k
Out[17]:
In [18]:
# 灰度变换
im = np.array(pil_im.copy().convert('L'))
Image.fromarray(255 - im) # 对图像反相处理
Out[18]:
In [19]:
Image.fromarray(uint8((100.0/255.0)*im + 100)) # 将像素值变换到 100..200
Out[19]:
In [20]:
Image.fromarray(uint8(255.0 * (im/255.0)**2)) # 将图像像素求平方后的图像
Out[20]:
In [21]:
def imresize(im,sz):
pil_im = Image.fromarray(uint8(im))
return array(pil_im.resize(sz))
In [22]:
def histeq(im,nbr_bins=256):
""" Histogram equalization of a grayscale image. """
if type(im) == Image.Image:
im = array(im)
# 计算直方图
imhist,bins = histogram(im.flatten(),nbr_bins,normed=True)
cdf = imhist.cumsum() # 计算分布函数
cdf = 255 * cdf / cdf[-1] # normalize 归一化
# 用积累分布函数的线性插值,计算新的像素值
im2 = interp(im.flatten(),bins[:-1],cdf)
im2 = Image.fromarray(uint8(im2.reshape(im.shape)))
return im2,cdf
In [23]:
# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
In [24]:
im = pil_im.copy()
im2, cdf = histeq(im)
In [25]:
axis('off')
title(u'原始图像', fontproperties=font)
imshow(im)
Out[25]:
In [26]:
axis('off')
title(u'直方图均衡化后的图像', fontproperties=font)
imshow(im2)
Out[26]:
In [27]:
axis('off')
title(u'原始直方图', fontproperties=font)
hist(array(im).flatten(), 128, normed=True)
show()
In [28]:
axis('off')
title(u'均衡化后的直方图', fontproperties=font)
hist(array(im2).flatten(), 128, normed=True)
show()
In [29]:
im = cv2.cvtColor(array(pil_im),cv2.COLOR_BGR2GRAY) # 转为灰度图像
imshow(cv2.equalizeHist(im))
Out[29]:
In [30]:
def compute_average(imlist):
""" Compute the average of a list of images. """
# open first image and make into array of type float
averageim = array(Image.open(imlist[0]), 'f')
skipped = 0
for imname in imlist[1:]:
try:
averageim += array(Image.open(imname))
except:
print imname + "...skipped"
skipped += 1
averageim /= (len(imlist) - skipped)
# return average as uint8
return array(averageim, 'uint8')
In [31]:
def pca(X):
""" 主成份分析
输入: X, 矩阵,每一行为一条训练数据(训练图像)
返回: 投影矩阵(按照维度的重要性排序),方差,均值
"""
# 获取维数
num_data,dim = X.shape
# 数据中心化
mean_X = X.mean(axis=0)
X = X - mean_X
if dim>num_data:
# PCA - 使用紧致技巧
M = dot(X,X.T) # 协方差矩阵 covariance matrix
e,EV = linalg.eigh(M) # 特征值和特征向量 eigenvalues and eigenvectors
tmp = dot(X.T,EV).T # this is the compact trick(紧致技巧)
V = tmp[::-1] # 由于最后的特征向量是我们需要的,所以需要进行逆转(reverse)
S = sqrt(e)[::-1] # 由于特征值是递增排序的,所以也要逆转
for i in range(V.shape[1]):
V[:,i] /= S
else:
# PCA - SVD 使用方法
U,S,V = linalg.svd(X)
V = V[:num_data] # 仅返回前 num_data 维的数据才合理
# 返回: 投影矩阵,方差,均值
return V,S,mean_X
In [32]:
import os
def get_imlist(path):
""" Returns a list of filenames for
all jpg images in a directory. """
return [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]
In [33]:
# 获取图像列表和他们的尺寸
imlist = get_imlist('a_thumbs')
im = array(Image.open(imlist[0])) # open one image to get the size
m, n = im.shape[:2] # get the size of the images
imnbr = len(imlist) # get the number of images
print "The number of images is %d" % imnbr
# Create matrix to store(存储) all flattened images(压平的图像)
immatrix = array([array(Image.open(imname)).flatten() for imname in imlist], 'f')
# PCA降维
V, S, immean = pca(immatrix)
# 显示图像(均值,和前 7 维图像
figure()
gray()
subplot(2, 4, 1)
axis('off')
imshow(immean.reshape(m, n))
for i in range(7):
subplot(2, 4, i+2)
imshow(V[i].reshape(m, n))
axis('off')
show()
主要函数说明:
cv2.PCACompute(data[, mean[, eigenvectors[, maxComponents ]]]) -> mean, eigenvectors
cv2.PCAComputeVar(data, retainedVariance[, mean[, eigenvectors ]]) -> mean, eigenvectors
cv2.PCAProject(data, mean, eigenvectors[, result ]) -> result
cv2.PCABackProject(data, mean, eigenvectors[, result ]) -> result
其中:
主要函数: cv2.PCACompute cv2.PCAComputeVar cv2.PCAProject cv2.PCABackProject
In [34]:
import cv2
mean, eig = cv2.PCACompute(immatrix,None)
# 显示图像(均值,和前 7 维图像
figure()
gray()
subplot(2, 4, 1)
axis('off')
imshow(mean.reshape(m, n))
for i in range(7):
subplot(2, 4, i+2)
imshow(eig[i].reshape(m, n))
axis('off')
show()
In [35]:
from PIL import Image
from pylab import *
from scipy.ndimage import filters
im = array(Image.open('empire.jpg'))
im2 = filters.gaussian_filter(im,5) # 标准差为5
axis('off')
imshow(im2)
Out[35]:
# 高斯模糊
cv2.GaussianBlur(src, ksize, sigmaX[, dst[, sigmaY[, borderType ]]]) -> dst
# 均值模糊
cv2.medianBlur(src, ksize[, dst ]) -> dst
# 通用模糊
cv2.filter2D(src, ddepth, kernel[, dst[, anchor[, delta[, borderType ]]]]) -> dst
# 获取 kernel
cv2.getGaussianKernel(ksize, sigma[, ktype ]) -> retval
cv2.getDerivKernels(dx, dy, ksize[, kx[, ky[, normalize[, ktype ]]]]) -> kx, ky
ksize : kernel size. ksize.width 和 ksize.height 可以不同,但必须为正的奇数,当其为0时从sigma计算得到
ddepth : 图像深度
通用模糊 filter2D:
$$ dst(x,y) = \sum_{ 0 \le x' < kernel.cols} \sum_{0 \le y' < kernel.rows} {kernel(x',y')*src(x + x' - anchor.x, y+y'-anchor.y)} $$主要函数: cv2.GaussianBlur cv2.medianBlur cv2.filter2D cv2.getGaussianKernel cv2.getDerivKernels
In [36]:
import cv2
im2 = cv2.GaussianBlur(im,(0,0),sigmaX=5,sigmaY=5)
axis('off')
imshow(im2)
Out[36]:
Prewitt 滤波器:
$$ D_x = \left[\begin{matrix}-1 & 0 & 1\\-1 & 0 & 1\\-1 & 0 & 1\end{matrix}\right] \\ D_y = \left[\begin{matrix}-1 & -1 & -1\\0 & 0 & 0\\1 & 1 & 1\end{matrix}\right] $$Sobel 滤波器:
$$ D_x = \left[\begin{matrix}-1 & 0 & 1\\-2 & 0 & 2\\-1 & 0 & 1\end{matrix}\right] \\ D_y = \left[\begin{matrix}-1 & -2 & -1\\0 & 0 & 0\\1 & 2 & 1\end{matrix}\right] $$
In [37]:
from scipy.ndimage import filters
from pylab import *
# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
im = array(Image.open('empire.jpg').convert('L'))
axis('off')
title(u'(a)原图', fontproperties=font)
imshow(im)
Out[37]:
In [38]:
# Sobel derivative filters
imx = zeros(im.shape)
filters.sobel(im, 1, imx)
axis('off')
title(u'(b)x方向差分', fontproperties=font)
imshow(imx)
Out[38]:
In [39]:
imy = zeros(im.shape)
filters.sobel(im, 0, imy)
axis('off')
title(u'(c)y方向差分', fontproperties=font)
imshow(imy)
Out[39]:
In [40]:
#mag = numpy.sqrt(imx**2 + imy**2)
mag = 255-numpy.sqrt(imx**2 + imy**2)
title(u'(d)梯度幅度', fontproperties=font)
axis('off')
imshow(mag)
Out[40]:
cv2.Sobel(src, ddepth, dx, dy[, dst[, ksize[, scale[, delta[, borderType ]]]]]) -> dst
cv2.Laplacian(src, ddepth[, dst[, ksize[, scale[, delta[, borderType ]]]]]) -> dst
ddepth -> 图像深度
主要函数: cv2.Sobel cv2.Scharr cv2.Laplacian
In [41]:
imx = cv2.Sobel(im,0,1,0)
axis('off')
imshow(imx)
Out[41]:
In [42]:
imy = cv2.Sobel(im,0,0,1)
axis('off')
imshow(imy)
Out[42]:
In [43]:
mag = cv2.Sobel(im,0,1,1)
axis('off')
imshow(mag)
Out[43]:
In [44]:
def imx(im, sigma):
imgx = zeros(im.shape)
filters.gaussian_filter(im, sigma, (0, 1), imgx)
return imgx
def imy(im, sigma):
imgy = zeros(im.shape)
filters.gaussian_filter(im, sigma, (1, 0), imgy)
return imgy
def mag(im, sigma):
# there's also gaussian_gradient_magnitude()
#mag = numpy.sqrt(imgx**2 + imgy**2)
imgmag = 255 - numpy.sqrt(imgx ** 2 + imgy ** 2)
return imgmag
im = array(Image.open('empire.jpg').convert('L'))
figure()
gray()
sigma = [2, 5, 10]
for i in sigma:
subplot(3, 4, 4*(sigma.index(i))+1)
axis('off')
imshow(im)
imgx=imx(im, i)
subplot(3, 4, 4*(sigma.index(i))+2)
axis('off')
imshow(imgx)
imgy=imy(im, i)
subplot(3, 4, 4*(sigma.index(i))+3)
axis('off')
imshow(imgy)
imgmag=mag(im, i)
subplot(3, 4, 4*(sigma.index(i))+4)
axis('off')
imshow(imgmag)
In [45]:
imx = cv2.GaussianBlur(im,(1,0),sigmaX=10,sigmaY=10)
axis('off')
imshow(imx)
Out[45]:
In [46]:
imy = cv2.GaussianBlur(im,(0,1),sigmaX=10,sigmaY=10)
axis('off')
imshow(imy)
Out[46]:
In [47]:
from scipy.ndimage import measurements, morphology
from pylab import *
im = array(Image.open('houses.png').convert('L'))
imshow(im)
axis('off')
title(u'原图', fontproperties=font)
Out[47]:
In [48]:
im = (im < 128)
labels, nbr_objects = measurements.label(im)
print "Number of objects:", nbr_objects
imshow(labels)
axis('off')
title(u'标记后的图', fontproperties=font)
Out[48]:
In [49]:
# morphology - opening to separate objects better
im_open = morphology.binary_opening(im, ones((9, 5)), iterations=2)
imshow(im_open)
axis('off')
title(u'开运算后的图像', fontproperties=font)
Out[49]:
In [50]:
labels_open, nbr_objects_open = measurements.label(im_open)
print "Number of objects:", nbr_objects_open
imshow(labels_open)
axis('off')
title(u'开运算后进行标记后的图像', fontproperties=font)
Out[50]:
cv2.erode cv2.dilate cv2.morphologyEx
# 腐蚀
cv2.erode(src, kernel[, dst[, anchor[, iterations[, borderType[, borderValue]]]]]) → dst
# 膨胀
cv2.dilate(src, kernel[, dst[, anchor[, iterations[, borderType[, borderValue]]]]]) → dst
# 形态学操作
cv2.morphologyEx(src, op, kernel[, dst[, anchor[, iterations[, borderType[, borderValue]]]]]) → dst
腐蚀:
$$ dst(x,y) = \max_{(x',y'):element(x',y') \ne 0}{src(x+x',y+y')} $$膨胀:
$$ dst(x,y) = \min_{(x',y'):element(x',y') \ne 0}{src(x+x',y+y')} $$从OpenCV 3.0 开始有函数:
int connectedComponents(InputArray image, OutputArray labels, int connectivity=8, int ltype=CV_32S)
int connectedComponentsWithStats(InputArray image, OutputArray labels, OutputArray stats, OutputArray centroids, int connectivity=8, int ltype=CV_32S)
现在也可以用 findContours 实现标记
此外常用算法有两遍扫描(Two-Pass)或平面区域填充(Seed-Filling),或满水填充算法(cv2.floodFill)
cv2.floodFill(image, mask, seedPoint, newVal[, loDiff[, upDiff[, flags]]]) → retval, rect)
a、这里的扫描指的是按行或按列访问以便图像的所有像素，本文算法采用的是按行扫描方式；
b、图像记为B，为二值图像：前景像素（pixel value = 1），背景像素（pixel value = 0）
c、label从2开始计数；
d、像素相邻关系：4-领域、8-领域，本文算法采用4-邻域；
（1）第一次扫描：
访问当前像素B(x,y)，如果B(x,y) == 1：
a、如果B(x,y)的领域中像素值都为0，则赋予B(x,y)一个新的label：
label += 1， B(x,y) = label；
b、如果B(x,y)的领域中有像素值 > 1的像素Neighbors：
1）将Neighbors中的最小值赋予给B(x,y):
B(x,y) = min{Neighbors}
2）记录Neighbors中各个值（label）之间的相等关系，即这些值（label）同属同一个连通区域；
labelSet[i] = { label_m, .., label_n }，labelSet[i]中的所有label都属于同一个连通区域
（注：这里可以有多种实现方式，只要能够记录这些具有相等关系的label之间的关系即可）
（2）第二次扫描：
访问当前像素B(x,y)，如果B(x,y) > 1：
a、找到与label = B(x,y)同属相等关系的一个最小label值，赋予给B(x,y)；
完成扫描后，图像中具有相同label值的像素就组成了同一个连通区域。
// Connected Component Analysis/Labeling By Two-Pass Algorithm
// Author: www.icvpr.com
// Blog : http://blog.csdn.net/icvpr
#include <iostream>
#include <string>
#include <list>
#include <vector>
#include <map>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
void icvprCcaByTwoPass(const cv::Mat& _binImg, cv::Mat& _lableImg)
{
// connected component analysis (4-component)
// use two-pass algorithm
// 1. first pass: label each foreground pixel with a label
// 2. second pass: visit each labeled pixel and merge neighbor labels
//
// foreground pixel: _binImg(x,y) = 1
// background pixel: _binImg(x,y) = 0
if (_binImg.empty() ||
_binImg.type() != CV_8UC1)
{
return ;
}
// 1. first pass
_lableImg.release() ;
_binImg.convertTo(_lableImg, CV_32SC1) ;
int label = 1 ; // start by 2
std::vector<int> labelSet ;
labelSet.push_back(0) ; // background: 0
labelSet.push_back(1) ; // foreground: 1
int rows = _binImg.rows - 1 ;
int cols = _binImg.cols - 1 ;
for (int i = 1; i < rows; i++)
{
int* data_preRow = _lableImg.ptr<int>(i-1) ;
int* data_curRow = _lableImg.ptr<int>(i) ;
for (int j = 1; j < cols; j++)
{
if (data_curRow[j] == 1)
{
std::vector<int> neighborLabels ;
neighborLabels.reserve(2) ;
int leftPixel = data_curRow[j-1] ;
int upPixel = data_preRow[j] ;
if ( leftPixel > 1)
{
neighborLabels.push_back(leftPixel) ;
}
if (upPixel > 1)
{
neighborLabels.push_back(upPixel) ;
}
if (neighborLabels.empty())
{
labelSet.push_back(++label) ; // assign to a new label
data_curRow[j] = label ;
labelSet[label] = label ;
}
else
{
std::sort(neighborLabels.begin(), neighborLabels.end()) ;
int smallestLabel = neighborLabels[0] ;
data_curRow[j] = smallestLabel ;
// save equivalence
for (size_t k = 1; k < neighborLabels.size(); k++)
{
int tempLabel = neighborLabels[k] ;
int& oldSmallestLabel = labelSet[tempLabel] ;
if (oldSmallestLabel > smallestLabel)
{
labelSet[oldSmallestLabel] = smallestLabel ;
oldSmallestLabel = smallestLabel ;
}
else if (oldSmallestLabel < smallestLabel)
{
labelSet[smallestLabel] = oldSmallestLabel ;
}
}
}
}
}
}
// update equivalent labels
// assigned with the smallest label in each equivalent label set
for (size_t i = 2; i < labelSet.size(); i++)
{
int curLabel = labelSet[i] ;
int preLabel = labelSet[curLabel] ;
while (preLabel != curLabel)
{
curLabel = preLabel ;
preLabel = labelSet[preLabel] ;
}
labelSet[i] = curLabel ;
}
// 2. second pass
for (int i = 0; i < rows; i++)
{
int* data = _lableImg.ptr<int>(i) ;
for (int j = 0; j < cols; j++)
{
int& pixelLabel = data[j] ;
pixelLabel = labelSet[pixelLabel] ;
}
}
}
In [51]:
def icvprCcaByTwoPass(im):
# 第一遍扫描
lableImg = np.array(im,dtype=np.float32) # 因为标签的值可能超过 255
rows, cols = im.shape
label = 1 # 从2开始
labelSet = []
labelSet.append(0) # 背景 0
labelSet.append(1) # 前景 1
for row in range(1,rows-1):
perRow = lableImg[row-1,:]
curRow = lableImg[row,:]
for col in range(1,cols-1):
if curRow[col] == 1:
neighborLabels = []
leftPixel = curRow[col-1]
upPixel = perRow[col]
if leftPixel > 1:
neighborLabels.append(leftPixel)
if upPixel > 1:
neighborLabels.append(upPixel)
if neighborLabels == []:
label += 1
labelSet.append(label) # 应用一个新的标签
curRow[col] = label
#labelSet[label] = label # 是不是多余? 不多余?
else:
neighborLabels = sort(neighborLabels)
smallestLabel = neighborLabels[0]
smallestLabel = int(smallestLabel)
curRow[col] = smallestLabel
# 保存等价的标签
for k in range(1,len(neighborLabels)):
tempLabel = neighborLabels[k]
tempLabel = int(tempLabel)
oldSmallestLabel = labelSet[tempLabel]
oldSmallestLabel = int(oldSmallestLabel)
if oldSmallestLabel > smallestLabel:
labelSet[oldSmallestLabel] = smallestLabel
labelSet[tempLabel] = smallestLabel
elif oldSmallestLabel < smallestLabel:
labelSet[smallestLabel] = oldSmallestLabel
# 更新等价的标签
# 确保 labelSet[Label] 中保存的是与 标签Label等价的最小Label
for i in range(2,len(labelSet)):
curLabel = labelSet[i]
preLabel = labelSet[curLabel]
while preLabel != curLabel:
curLabel = preLabel
preLabel = labelSet[preLabel]
labelSet[i] = curLabel
# 第二遍扫描更新图像
for row in range(rows):
for col in range(cols):
lableImg[row,col] = labelSet[int(lableImg[row,col])]
# 或者如下
# lableImg = np.array([labelSet[int(x)] for x in lableImg])
return lableImg, len(set(labelSet)) - 2 # 因为 0,1 两个不是标签,因此标签数量应该是总标签数量减去2
In [52]:
im = array(Image.open('houses.png').convert('L'))
im = cv2.threshold(im,128,1,type=cv2.THRESH_BINARY_INV)[1]
lableImg, labelNum = icvprCcaByTwoPass(im)
print(labelNum)
imshow(lableImg)
Out[52]:
// Connected Component Analysis/Labeling By Seed-Filling Algorithm
// Author: www.icvpr.com
// Blog : http://blog.csdn.net/icvpr
#include <iostream>
#include <string>
#include <list>
#include <vector>
#include <map>
#include <stack>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
void icvprCcaBySeedFill(const cv::Mat& _binImg, cv::Mat& _lableImg)
{
// connected component analysis (4-component)
// use seed filling algorithm
// 1. begin with a foreground pixel and push its foreground neighbors into a stack;
// 2. pop the top pixel on the stack and label it with the same label until the stack is empty
//
// foreground pixel: _binImg(x,y) = 1
// background pixel: _binImg(x,y) = 0
if (_binImg.empty() ||
_binImg.type() != CV_8UC1)
{
return ;
}
_lableImg.release() ;
_binImg.convertTo(_lableImg, CV_32SC1) ;
int label = 1 ; // start by 2
int rows = _binImg.rows - 1 ;
int cols = _binImg.cols - 1 ;
for (int i = 1; i < rows-1; i++)
{
int* data= _lableImg.ptr<int>(i) ;
for (int j = 1; j < cols-1; j++)
{
if (data[j] == 1)
{
std::stack<std::pair<int,int>> neighborPixels ;
neighborPixels.push(std::pair<int,int>(i,j)) ; // pixel position: <i,j>
++label ; // begin with a new label
while (!neighborPixels.empty())
{
// get the top pixel on the stack and label it with the same label
std::pair<int,int> curPixel = neighborPixels.top() ;
int curX = curPixel.first ;
int curY = curPixel.second ;
_lableImg.at<int>(curX, curY) = label ;
// pop the top pixel
neighborPixels.pop() ;
// push the 4-neighbors (foreground pixels)
if (_lableImg.at<int>(curX, curY-1) == 1)
{// left pixel
neighborPixels.push(std::pair<int,int>(curX, curY-1)) ;
}
if (_lableImg.at<int>(curX, curY+1) == 1)
{// right pixel
neighborPixels.push(std::pair<int,int>(curX, curY+1)) ;
}
if (_lableImg.at<int>(curX-1, curY) == 1)
{// up pixel
neighborPixels.push(std::pair<int,int>(curX-1, curY)) ;
}
if (_lableImg.at<int>(curX+1, curY) == 1)
{// down pixel
neighborPixels.push(std::pair<int,int>(curX+1, curY)) ;
}
}
}
}
}
}
In [53]:
def icvprCcaBySeedFill(img):
"""
使用如下算法:
1. begin with a foreground pixel and push its foreground neighbors into a stack;
2. pop the top pixel on the stack and label it with the same label until the stack is empty
foreground pixel: img(x,y) = 1
background pixel: img(x,y) = 0
"""
lableImg = np.array(im,dtype=np.float32) # 因为标签的值可能超过 255
label = 1 # 标签从 2 开始
rows, cols = img.shape
for row in range(1,rows-1):
data = lableImg[row,:]
for col in range(1,cols-1):
if data[col] == 1:
neighborPixels = []
neighborPixels.append((row,col)) # 起始像素位置
label += 1 # 开始一个新的标签
while not neighborPixels == []:
# 获取堆栈顶层的像素位置并将其对应的像素标记为同一个 label
curPixel = neighborPixels[-1]
curRow, curCol = curPixel
lableImg[curPixel] = label
# 弹出堆栈顶层的像素位置
neighborPixels.pop()
# 将前景(像素值为1)中的 4-邻接 像素的位置压入堆栈
neighbors_4 = []
if curCol-1 > -1:
neighbors_4.append((curRow,curCol-1)) # 左邻接像素
if curCol+1 < cols:
neighbors_4.append((curRow,curCol+1)) # 右邻接像素
if curRow-1 > -1:
neighbors_4.append((curRow-1,curCol)) # 上邻接像素
if curRow+1 < rows:
neighbors_4.append((curRow+1,curCol)) # 下邻接像素
for pos in neighbors_4:
if lableImg[pos] == 1:
neighborPixels.append(pos)
return lableImg, label-1 # 因为 1 不是标签所以标签数量应该是标签值减去1
In [54]:
im = array(Image.open('houses.png').convert('L'))
im = cv2.threshold(im,128,1,type=cv2.THRESH_BINARY_INV)[1]
lableImg, labelNum = icvprCcaBySeedFill(im)
print(labelNum)
imshow(lableImg)
Out[54]:
In [55]:
def denoise(im,U_init,tolerance=0.1,tau=0.125,tv_weight=100):
""" An implementation of the Rudin-Osher-Fatemi (ROF) denoising model
using the numerical procedure presented in Eq. (11) of A. Chambolle
(2005). Implemented using periodic boundary conditions.
Input: noisy input image (grayscale,灰度), initial guess for U, weight of
the TV-regularizing term(正则项权值), steplength, tolerance for the stop criterion
Output: denoised(去噪) and detextured image, texture residual(残留). """
m,n = im.shape # 噪声图像大小
# 初始化
U = U_init
Px = zeros((m, n)) # 对偶域的x分量
Py = zeros((m, n)) # 对偶域的y分量
error = 1
while (error > tolerance):
Uold = U
# 原始变量的梯度
GradUx = roll(U,-1,axis=1)-U # 变量U梯度的 x 分量
GradUy = roll(U,-1,axis=0)-U # 变量U梯度的 y 分量
# 更新对偶变量
PxNew = Px + (tau/tv_weight)*GradUx # non-normalized update of x-component (dual 对偶域)
PyNew = Py + (tau/tv_weight)*GradUy # non-normalized update of y-component (dual)
NormNew = maximum(1,sqrt(PxNew**2+PyNew**2))
Px = PxNew/NormNew # update of x-component (dual)
Py = PyNew/NormNew # update of y-component (dual)
# update the primal variable(原始变量)
RxPx = roll(Px,1,axis=1) # right x-translation of x-component
RyPy = roll(Py,1,axis=0) # right y-translation of y-component
DivP = (Px-RxPx)+(Py-RyPy) # divergence(散度) of the dual field.
U = im + tv_weight*DivP # update of the primal variable
# update of error
error = linalg.norm(U-Uold)/sqrt(n*m);
return U,im-U # denoised image and texture residual
In [56]:
from numpy import random
# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
# create synthetic image with noise
im = zeros((500,500))
im[100:400,100:400] = 128
im[200:300,200:300] = 255
im = im + 30*random.standard_normal((500,500))
U,T = denoise(im,im)
G = filters.gaussian_filter(im,10)
# save the result
#imsave('synth_original.pdf',im)
#imsave('synth_rof.pdf',U)
#imsave('synth_gaussian.pdf',G)
# plot
figure()
gray()
subplot(1,3,1)
imshow(im)
#axis('equal')
axis('off')
title(u'原噪声图像', fontproperties=font)
subplot(1,3,2)
imshow(G)
#axis('equal')
axis('off')
title(u'高斯模糊后的图像', fontproperties=font)
subplot(1,3,3)
imshow(U)
#axis('equal')
axis('off')
title(u'ROF降噪后的图像', fontproperties=font)
Out[56]:
In [ ]: