图像处理 高斯拟合

2017-06-21 09:25:56 wangkun1340378 阅读数 5614

目的: matlab 图像拟合

效果图:

 


利用函数: surf

代码:

theat=1:1:300;
fai=1:1:400;
[t,f]=meshgrid(theat,fai);
z=imread('test.jpg');
z=imresize(z,[300 400]);   %%大小resize成300X400
z=z(:,:,1);
z=z';

surf(t,f,z,'edgecolor','none')

alpha(1)     %%这个1代表透明度 可选范围是0-1 1是完全不透明 




2018-06-13 12:56:20 ma7856728 阅读数 8810

在边缘检测中总会提取出不连续点,或伪轮廓。在这种情况下需要拟合出目标的轮廓,这样可以找到轮廓的数学表达式为后续的特征选取打下基础。博主用coins图像为例,用椭圆方程进行拟合,做出如下实验。

1、原图二值化

2、边缘检测(sobel算子)

3、填补孔洞

4、标记连通域

5、找到每个连通域坐标

6、用每个连通域坐标拟合出椭圆方程

7、在二值图像中画出每个椭圆函数

%%图像边缘检测和拟合轮廓
clc
clear
close all
%% 读取图像
I = imread('coins.png');
I = im2bw(I);                                   %二值化
I = imfill(I,'holes');                          %填补孔洞
[M,N] = size(I);
figure(1),imshow(I);title('原图');hold on
%% 选取待拟合坐标
% conicP = ginput(15);            
bw1 = edge(I,'sobel');                           %边缘检测
figure(2),imshow(bw1);title('边缘检测');
[L,num] = bwlabel(bw1);                           %标签
for i = 1:num
    [row,col] = find(L == i);
    conicP = zeros(length(row),2);
    conicP(:,1) = col;
    conicP(:,2) = row;
    figure(1),plot(conicP(:,1)', conicP(:,2)', 'xr');     %drawing sample points
%% 自定义椭圆函数拟合
    a0 = [1 1 1 1 1 1];
    f = @(a,x)a(1)*x(:,1).^2+a(2)*x(:,2).^2+a(3)*x(:,1).*x(:,2)+a(4)*x(:,1)+a(5)*x(:,2)+a(6);%建立方程
    p = nlinfit(conicP , zeros(size(conicP, 1), 1), f,[1 2 3 4 5 6]);
    syms x y
    conic = p(1)*x^2+p(2)*y^2+p(3)*x*y+p(4)*x+p(5)*y+p(6);
%% 在原图上显示拟合结果
    c = ezplot(conic,[0,N],[0,M]);
    figure(1),set(c, 'Color', 'Blue','LineWidth',2);
   
end
 hold off

实验结果图像如下,拟合结果为蓝线,原坐标为红点


2018-12-21 15:44:23 yql_617540298 阅读数 1044

一、最小二乘法拟合曲线

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

#自定义函数 e指数形式
def func(x, a, b,c):
    return a*np.sqrt(x)*(b*np.square(x)+c)

#定义x、y散点坐标
x = [10,20,30,40,50,60,70,80]
x = np.array(x)
y = [158,455,265,152,263,813,562,126]
y = np.array(y)

#非线性最小二乘法拟合
popt, pcov = curve_fit(func, x, y)
#获取popt里面是拟合系数
print(popt)
a = popt[0]
b = popt[1]
c = popt[2]
yvals = func(x,a,b,c) #拟合y值
print('popt:', popt)
print('系数a:', a)
print('系数b:', b)
print('系数c:', c)
print('系数pcov:', pcov)
print('系数yvals:', yvals)
#绘图
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()

二、高斯分布拟合曲线

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


#自定义函数 e指数形式
def func(x, a,u, sig):
    return  a*(np.exp(-(x - u) ** 2 /(2* sig **2))/(math.sqrt(2*math.pi)*sig))*(431+(4750/x))


#定义x、y散点坐标
x = [10,20,30,40,50,60,70,80]
x=np.array(x)
# x = np.array(range(20))
print('x is :\n',x)
y = [158,455,265,152,263,813,562,126]
y = np.array(y)
print('y is :\n',y)

popt, pcov = curve_fit(func, x, y,p0=[3.1,4.2,3.3])
#获取popt里面是拟合系数
a = popt[0]
u = popt[1]
sig = popt[2]


yvals = func(x,a,u,sig) #拟合y值
print(u'系数a:', a)
print(u'系数u:', u)
print(u'系数sig:', sig)

#绘图
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()

三、多项式拟合曲线

import numpy as np
import matplotlib.pyplot as plt

#定义x、y散点坐标
#x = [10,20,30,40,50,60,70,80]
#y = [158,455,265,152,263,813,562,126]
#x = [16,32,48,64,80,96,112,128,144,160,176,192,208,224,240,256,272,288,304]
x = [0,16,32,48,64,80,96,112,128,144,160,176,192,208,224,240,256,272,288,304,320,336,352,368,384,400,416,432,448,464,480,496,]
x = np.array(x)
#y = [506,506,506,506,506,506,506,506,506,505,505,505,505,505,505,505,506,505,504]
y = [340,338,345,348,348,349,50,350,350,350,350,350,350,350,350,350,350,349,348,348,348,347,347,348,348,347,347,347,347,348,347,346]
y = np.array(y)
#用3次多项式拟合
f1 = np.polyfit(x, y, 3)
p1 = np.poly1d(f1)
yvals = p1(x)#拟合y值

#x1 = [32,48,64,80,96,112,128,144,160,176,192,208,224,240,256,272,288,304,320]
x1 = [16,32,48,64,80,96,112,128,144,160,176,192,208,224,240,256,272,288,304,320,336,352,368,384,400,416,432,448,464,480,496,512]
x1 = np.array(x1)
#y1 = [529,528,528,529,529,529,529,529,529,529,528,528,528,528,528,528,528,528,529]
y1 = [370,367,376,378,378,379,379,380,380,380,379,379,379,379,378,378,378,378,377,376,376,375,375,372,372,372,372,375,375,372,372,373]
y1 = np.array(y1)
f2 = np.polyfit(x1,y1,3)
p2 = np.poly1d(f2)
yvals_2 = p2(x1)
#也可使用yvals=np.polyval(f1, x)
#绘图
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')

plot3 = plt.plot(x1, y1, 's', label = 'original values')
plot4 = plt.plot(x1, yvals_2, 'r',label='',color='green')

#plt.axis('off')
#plt.legend(loc=4) #指定legend的位置右下角
#plt.title('polyfitting')
#plt.show()
plt.savefig("F:/a.jpg")

2019-11-25 09:58:01 m0_37816922 阅读数 321

通过R语言处理高斯光束

图片加载与显示

R语言中也有不少图像处理包,著名的magick就提供了R语言的接口。但是magick包更像是一个代码版的PS,可以实现诸多高级功能,但过于完整的代码封装使得一些基础操作反而难以施展。所以我们使用imager包。
本实验在RStudio中运行。
首先,安装并导入包。
然后,通过load.image读取图片赋值给变量,在R语言中推荐使用<-进行赋值,赋值之后可以在右侧Environment选项卡中查看数据区,可以看到我们赋值的img是一个Large cimg类型。
通过dim函数可以查看数据的维度,可以看到该数据类型有4个维度,前两个分别代表图片的长、宽;第三个代表帧数;第四个代表颜色通道,其值为3,表示这是一个rgb图片。
读取图片之后通过plot命令可以方便地查看图片数据。

> install.packages(imager)
> library(imager)
> img <- load.image("1.bmp")
> dim(img)
[1] 640 480   1   3
> plot(img)

在使用plot之后,可以在右侧的Plots选项卡中查看刚刚画出的数据,如下图所示。
在这里插入图片描述

图像截取

为了便于处理,一般需要将图片转换为矩阵格式,而对于rgb图像而言,还需要将通道统一,即先通过grayscale转成灰度图像,然后通过as.matrix转成矩阵。
通过image函数可以查看矩阵的数据,其输入参数为x轴坐标取值、y轴坐标取值。

> gray <- grayscale(img)
> imgMat = as.matrix(gray)
> dim(imgMat)
[1] 640 480
> m = dim(imgMat)[1]
> n = dim(imgMat)[2]
> x <- 1:m
> y <- 1:n
> image(x,y,imgMat)

在这里插入图片描述

我们可以看到矩阵中大量的数据是无效的背景,需要对矩阵进行截取。在得到矩阵的image图像之后,可以通过locator命令交互式选取坐标位置,从而可以精确地截取矩阵。其方法为,输入命令之后,将鼠标放在图像上,当鼠标光标变为十字时点击,R会记录下鼠标点击的位置。
然后通过[:,:]的索引方法,对矩阵进行截取。

> unlist(locator(2))
      x1       x2       y1       y2 
356.8567 422.5432 186.0054 246.7325
> roi = imgMat[356:422,186:246]
> image(roi)		#当不输入坐标时,其横纵坐标默认归一化

在这里插入图片描述

显示强度

伪彩图虽然也能反应出光斑的强度信息,但总体来说不够直观,而是依赖于人的主观感觉。在R语言中,graphics包提供了3d图像绘制函数persp,其输入参数分别为x、y、z,所以我们需要将矩阵拆分成三个向量,分别代表x轴、y轴以及图片的像素点强度。

> x <- 1:dim(roi)[1]
> y <- 1:dim(roi)[2]
> persp(x,y,z)
> persp(x,y,roi,theta=30,phi=50,axes=TRUE,shade=0.2,col="#00AAFF")

其中,theta,phi这两个参数的作用是调整视角,shade为阴影,col为颜色,其输入的值为16进制的rgb。

在这里插入图片描述

高斯拟合

由于矩阵是二维数据,所以我们首先通过每列选取最大值的方式得到一个数组。R语言中不直接提供计算矩阵每一行最大值的函数,但是提供了apply,可以指定矩阵的行操作或者列操作。其输入变量位待处理矩阵,行列标识以及处理函数。行列标识中,1代表行处理,2代表列处理。

> arr = apply(roi,1,max)	#对每一行执行max操作,max即选取最大值
> plot(arr)					#画图

在这里插入图片描述

最后,就是最重要的数据拟合。在R语言中,通过nlm函数进行非线性拟合,其输入参数分别为拟合函数表达式以及拟合参数的初始值。其中,拟合函数表达式中的变量需要预先定义好。
对于高斯函数,其表达式为

y=aexp((xbc)2)y=a⋅\exp{(−(\frac{x−b}{c})^2​)}

其中,y即我们之前求取最大值得到的arr,x则为与arr等长的自然数列。a的值表示该函数的最大值;b表示其中心值,c表示当y值降到1e\frac{1}{e}分之一处时x距离中心的位置。初始值可根据此选取。

> x <- 1:length(arr)		#建立x
> a0 <- max(arr)			#a的初始值
> b0 <- length(arr)/2	#b的初始值,根据图像可知,b的初始值在数据的中心位置
> c0 <- length(arr)/4	#c的初始值
> fit <- nlm(arr~~a*exp(-((x-b)/c)^2),start=list(a=a0,b=b0,c=c0))
> fit
Nonlinear regression model
  model: arr ~ a * exp(-((x - b)/c)^2)
   data: parent.frame()
      a       b       c 
 0.2285 35.7438 19.7131 				#此即为拟合参数
 residual sum-of-squares: 0.005557

Number of iterations to convergence: 5 
Achieved convergence tolerance: 5.998e-06

通过predict函数可以提取fit中的模型,并根据输入的x得出拟合之后的y值。我们可以据此画出原始数据与拟合函数的图像。

> y = predict(fit,list(x=x))
> plot(x,arr)
> lines(x,y,col='red')

在这里插入图片描述

2018-03-19 15:47:39 Sirius_007 阅读数 8049

1.最小二乘拟合

    最简单,但是受离群点的影响比较大,鲁棒性不强。

2.梯度下降法

3.高斯牛顿、列-马算法