2017-02-09 11:24:47 lpsl1882 阅读数 1217

空间域上,可以和频域一样进行卷积逆滤波操作。其方法是展开图像为一列,构建卷积模板矩阵,这样卷积操作就变成了矩阵乘法。我们可以用最小二乘法来,已知卷积图像和卷积模板来求出原始图像。空间域最小二乘逆滤波是病态问题,缺点是卷积矩阵非常稀疏和巨大,模非常小,一般需要进行约束。
我们现在有一个均值卷积模板blurKernel,对图像X进行卷积,生成一个模糊图像B,可以写成blurKernelX+N=B,其中*是卷积操作,N是加性噪声。展开后我们得到卷积矩阵A,以及B图像的展开Bravel,那么有

AXravel+Nravel=BravelXravel=(ATA+λI)1ATBravel

假设B是4x4矩阵,比如15913261014371115481216展开后为Bravel=[1,2,3,4,516]T;假设A的卷积是3x3的均值模板,不考虑边界,A是6x16矩阵,那么A=1/9......1/91/901/91/91/901/91/91/90000
可以想象到,这些矩阵非常稀疏和巨大。我们可以用稀疏矩阵库来进行计算,其中进行矩阵求逆是最为耗时耗力的一步。

import numpy as np
import os, string
from matplotlib import pyplot as plt
import scipy as sp
import cv2
img = cv2.imread('camera.jpg')
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
img = cv2.resize(img,(0,0),fx=0.5,fy=0.5)
print img.shape
from __future__ import division
blurKernel = np.ones((3,3))/9

from scipy import signal#warning
blurImg = signal.convolve2d(img, blurKernel, 'same','symm')

import itertools
from scipy.sparse import csc_matrix,lil_matrix
from scipy.sparse import linalg as sppl
import numpy.linalg
A=lil_matrix(((img.shape[0]-3)*(img.shape[1]-3),img.shape[0]*img.shape[1]))
print A.shape
ind = -1
for j in xrange(0,np.int32(img.shape[1])-3):
    for i in xrange(0,np.int32(img.shape[0])-3):
        ind +=1
        #index=[]
        ravel = []
        for e in itertools.product(range(j,j+3),range(i,i+3)):
            #index.append(e)
            ravel.append(e[1]+e[0]*img.shape[0])
        A[ind,[ravel]]=1./9.
 A=csc_matrix(A) 

from scipy import sparse
reg = 0.001
X = A.T.dot(A)
X=X+reg*sparse.eye(X.shape[0])
X=csc_matrix.dot(csc_matrix(Xinv), A.T)
padBlur = csc_matrix(blurImg[0:np.int32(img.shape[0])-3,0:np.int32(img.shape[1])-3].ravel()).T
print X.shape,padBlur.shape
X=csc_matrix.dot(X,padBlur)
res=X.todense()
res = res.reshape(blurImg.shape)
fig = plt.figure()
ax = fig.add_subplot(131)  
ax.imshow(img,cmap='gray')
ax = fig.add_subplot(132)  
ax.imshow(blurImg,cmap='gray')
ax = fig.add_subplot(133)  
ax.imshow(res,cmap='gray')
plt.show()

计算花费半个小时。相比频域逆滤波,时间和存储空间代价都很大,一般不推荐直接的空间域逆滤波。可以采用迭代优化的方法。
这里写图片描述
Fig. (a)原图 (b)均值滤波图像 (c)逆滤波图像

2017-02-09 19:51:25 lpsl1882 阅读数 1319

上一篇文章,我的空间域最小二乘逆滤波的时间、空间复杂度都非常高。其中求逆矩阵是消耗巨大的一步,这里用迭代优化解法展示了如何不用求逆矩阵来求解最小二乘逆滤波。
首先卷积图像的生成表示为AXravel+Nravel=Bravel,我们将问题转化为一个简单的带正则化的优化问题

J(Xravel)=argmin12(AXravelBravel)2+λ2|Xravel|2

这是一个无约束问题,我们直接用最速下降法求解:
JX=(ATA+λI)XravelATBravelXk+1ravel=Xkμ[(ATA+λI)XkravelATBravel]

其中X的初始值是一个随机矩阵。构建卷积矩阵A最耗时,而最速下降法就快多了。时间和空间复杂度都比直接最小二乘法要好的多。
下面是μ=0.7,λ=0.01epoch=200的结果
这里写图片描述
Fig (a)原图 (b)5x5均值滤波图像 (c)优化求解逆滤波图像

更新:
然而直接展开卷积作为矩阵乘法,卷积矩阵非常稀疏,即使使用第三方稀疏矩阵库,依然很麻烦。让我们解决这个问题。
解决方案有两个:

  • 加速卷积计算
  • 将转置卷积矩阵乘法变为普通卷积运算

第一个问题,在opencv和scipy中卷积已经获得了非常快的解法。opencv中加速卷积计算的要点是利用了计算机原理中的局部性原理。对于一个3x3的卷积,我们可以使用3个指针,第i个指针指向卷积区域的第i行,这样就形成了一个线扫描区域,每个指针顺序滑动指向下一个邻近数据地址,这样可以极大地利用计算机缓存,加快速度。

第二个问题,可以从CNN卷积神经网络的反向传播中获取思想。神经网络的反向传播,本质上是求解梯度下降法。而卷积层的反向梯度,即对卷积的求导操作,依然是一个卷积操作,其公式是

J=conv(σ,rot180(kernel))
也就是旋转了卷积模板,跟自相关计算类似。推导过程见:Convolutional Neural Networks backpropagation类比本文第一个公式,AX其实就是对X用K卷积模板(A是K卷积的展开)做卷积,ATΔX则是对AX矩阵乘法的求导,也就是XK卷积操作的求导,因此可以还原为卷积操作。
综合上述,我们得到的最优化超快速解法为:
Xk+1ravel=Xkμ{conv[conv(X,K)B,rot180(K)]+λX}

以下第一幅图是模糊图像,第二幅图是μ=1.0,λ=0.01epoch=200的逆滤波结果,第三幅图是原图
这里写图片描述

PS:这个算法相当于把一个随机图像变成卷积图像的逆滤波图像,与神经网络图像风格化neural style有共通之妙。

2015-08-10 14:30:47 fate08301017 阅读数 621

最近一直在进行图像处理方面的学习,在正式进行图像处理之前我们要学习许多预备的数学

知识,其中最小二乘法就是其中很重要的一种知识。

首先我们了解一下最小二乘法的原理与概念,最小二乘法(又称最小平方法)是一种数学优化技术。

它通过最小化误差的平方和寻找数据的最佳函数匹配。

利用最小二乘法可以进行最小二乘回归

       有一组数据(x1,y1,(x2,y2)......(xn,yn)

        我们假设这组数据拟合为y=ax这条直线,则用最小二乘法:

               

                           a应该取为满足以上公式取值最小时的值,

                           

这样我们就可以得到一条与这些已知数据拟合的比较好的直线y=ax

如果是一种非线性回归,比如y=a0+a1*x+a2*x^2

也可以看成线性的来描述,令x^2=x2,x=x1,1=x0,可以看成是一种对参数空间的线性,这样

就化成了y=a0*x0+a1*x1+a2*x2,这是一个多维的线性回归

如果有n维,我们可以有函数y=a1*x1+a2*x2+......+am*xm,我们可以

使用这样三个矩阵:

    A=  为参数矩阵    和   Y= 因变量观测值矩阵

以及X=自变量观测值矩阵

我们要求得最佳的参数即采用最小二乘法

误差为,我们取最小误差而得到参数值


我们如何由上式求出参数A呢,这里我们可以对误差函数求关于参数A的偏导

误差函数  可以化为

eA=求关于A的偏导可得

=0

由此可得A=

由此参数A可拟合出符合数据集的函数。

2015-06-18 17:00:56 hjimce 阅读数 13414

基于移动最小二乘的图像变形

原文地址:http://blog.csdn.net/hjimce/article/details/46550001

作者:hjimce

一、背景意义

写这篇博文是应为目前为止我看到了好多领域里的经典paper算法都有涉及到移动最小二乘(MLS)。可见这个算法非常重要,先来看一下它的相关经典应用:

1、图像变形。在图像处理领域paper:《Image Deformation Using Moving Least Squares》利用移动最小二乘的原理实现了图像的相关变形,而且这篇paper的引用率非常高,可以说是图像变形算法的经典算法,Siggraph上面的paper。


利用移动最小二乘实现图像变形

2、点云滤波。利用MLS实现点云滤波,是三维图像学点云处理领域的一大应用,我所知道点云滤波经典算法包括:双边滤波、MLS、WLOP。

3、Mesh Deformation。用这个算法实现三角网格模型的变形应用也是非常不错的,相关的paper《3D Deformation Using Moving Least Squares

OK,接着我就以《Image Deformation Using Moving Least Squares》算法为例,进行讲解基于移动最小二乘的图像变形算法实现。

二、算法实现

在这里我没有打算将算法原理的推导过程,直接讲算法的实现步骤公式。

这篇paper根据变换矩阵的不同,可以分为三种变形方法,分别是仿射变换、相似变换、刚性变换。其中刚性变换的效果是最好的,我这边从简单的讲,只讲仿射变换的变形算法实现:

问题:原图像的各个控制顶点坐标p,原图像上的像素点v的坐标。变形后图像的控制顶点位置q,求v在变形后图像中对应位置f(v)。

总计算公式为:


上面中lv(x)和f(v)是同一个函数。因为x就是我们输入的原图像的像素点坐标v。

因此我们的目标就是要知道p*,q*,变换矩阵M。这样输入一个参数x,我们就可以计算出它在变换后图像中的位置了。

OK,只要我们知道上面公式中,各个参数的计算方法,我们就可以计算出变形后图像对应的坐标点f(v)了。

1、权重w的计算方法为:


也就是计算v到控制顶点pi的距离倒数作为权重,参数a一般取值为1。

这一步实现代码如下:

	//计算各个控制顶点的权重,也就是计算点t到各个顶点的距离1/sqr(d)
	while(iter!=p.end())
	{
		double temp;
		if(iter->x!=t.x || iter->y!=t.y)
		    temp=1/((iter->x-t.x)*(iter->x-t.x)+(iter->y-t.y)*(iter->y-t.y));
		else//如果t为控制顶点,那么需要把该控制顶点的权重设置为无穷大
			temp=MAXNUM;
		w.push_back(temp);
		iter++;
	}
2、q*,p*的计算公式如下:

也就是计算控制顶点pi和qi的加权求和重心位置。

	double px=0,py=0,qx=0,qy=0,tw=0;
	while(iterw!=w.end())
	{
	px+=(*iterw)*(iter->x);//所有控制顶点p的加权位置
	py+=(*iterw)*(iter->y);
	qx+=(*iterw)*(iterq->x);//所有控制顶点q的加权位置
	qy+=(*iterw)*(iterq->y);
	tw+=*iterw;//总权重
	iter++;
    iterw++;
	iterq++;
	}
	pc.x=px/tw;
	pc.y=py/tw;
	qc.x=qx/tw;
	qc.y=qy/tw;
3、仿射变换矩阵M的计算公式如下:



只要把相关的参数都带进去就可以计算了。

最后贴一些完整的MLS源代码:

//输入原图像的t点,输出变形后图像的映射点f(v)
MyPoint CMLSDlg::MLS(const MyPoint& t)
{
	if(p.empty())//原图像的控制顶点p,与输入点t为同一副图像坐标系下
		return t;
	MyPoint fv;
	double A[2][2],B[2][2],M[2][2];
	iter=p.begin();
	w.erase(w.begin(),w.end());
	//计算各个控制顶点的权重,也就是计算点t到各个顶点的距离1/sqr(d)
	while(iter!=p.end())
	{
		double temp;
		if(iter->x!=t.x || iter->y!=t.y)
		    temp=1/((iter->x-t.x)*(iter->x-t.x)+(iter->y-t.y)*(iter->y-t.y));
		else//如果t为控制顶点,那么需要把该控制顶点的权重设置为无穷大
			temp=MAXNUM;
		w.push_back(temp);
		iter++;
	}
	vector<double>::iterator iterw=w.begin();
	vector<MyPoint>::iterator iterq=q.begin();//q为目标图像的控制点的位置,我们的目标是找到t在q中的对应位置
	iter=p.begin();
	MyPoint pc,qc;
	double px=0,py=0,qx=0,qy=0,tw=0;
	while(iterw!=w.end())
	{
	px+=(*iterw)*(iter->x);//所有控制顶点p的加权位置
	py+=(*iterw)*(iter->y);
	qx+=(*iterw)*(iterq->x);//所有控制顶点q的加权位置
	qy+=(*iterw)*(iterq->y);
	tw+=*iterw;//总权重
	iter++;
    iterw++;
	iterq++;
	}
	pc.x=px/tw;
	pc.y=py/tw;
	qc.x=qx/tw;
	qc.y=qy/tw;
	iter=p.begin();
	iterw=w.begin();
	iterq=q.begin();
	for(int i=0;i<2;i++)
    for(int j=0;j<2;j++)
	{
	A[i][j]=0;
	B[i][j]=0;
	M[i][j]=0;
	}
	while(iter!=p.end())
	{

	double P[2]={iter->x-pc.x,iter->y-pc.y};
	double PT[2][1];
	PT[0][0]=iter->x-pc.x;
	PT[1][0]=iter->y-pc.y;
	double Q[2]={iterq->x-qc.x,iterq->y-qc.y};
	double T[2][2];
   
    T[0][0]=PT[0][0]*P[0];
    T[0][1]=PT[0][0]*P[1];
	T[1][0]=PT[1][0]*P[0];
    T[1][1]=PT[1][0]*P[1];

	for(int i=0;i<2;i++)
    for(int j=0;j<2;j++)
	{
	A[i][j]+=(*iterw)*T[i][j];
	}
	T[0][0]=PT[0][0]*Q[0];
    T[0][1]=PT[0][0]*Q[1];
	T[1][0]=PT[1][0]*Q[0];
    T[1][1]=PT[1][0]*Q[1];
    
	for(int i=0;i<2;i++)
    for(int j=0;j<2;j++)
	{
	B[i][j]+=(*iterw)*T[i][j];
	}

	iter++;
    iterw++;
	iterq++;
	}
	//cvInvert(A,M);
    double det=A[0][0]*A[1][1]-A[0][1]*A[1][0];
	if(det<0.0000001)
	{
		fv.x=t.x+qc.x-pc.x;
		fv.y=t.y+qc.y-pc.y;
		return fv;
	}
    double temp1,temp2,temp3,temp4;
	temp1=A[1][1]/det;
	temp2=-A[0][1]/det;
	temp3=-A[1][0]/det;
	temp4=A[0][0]/det;
	A[0][0]=temp1;
    A[0][1]=temp2;
	A[1][0]=temp3;
	A[1][1]=temp4;


	M[0][0]=A[0][0]*B[0][0]+A[0][1]*B[1][0];
	M[0][1]=A[0][0]*B[0][1]+A[0][1]*B[1][1];
	M[1][0]=A[1][0]*B[0][0]+A[1][1]*B[1][0];
	M[1][1]=A[1][0]*B[0][1]+A[1][1]*B[1][1];

	double V[2]={t.x-pc.x,t.y-pc.y};
	double R[2][1];
  
	R[0][0]=V[0]*M[0][0]+V[1]*M[1][0];//lv(x)总计算公式
	R[1][0]=V[0]*M[0][1]+V[1]*M[1][1];
	fv.x=R[0][0]+qc.x;
	fv.y=R[1][0]+qc.y;

	return fv;
}


调用方法示例:

    int i=0,j=0;
	dImage=cvCreateImage(cvSize(2*pImage->width,2*pImage->height),pImage->depth,pImage->nChannels);//创建新的变形图像
	cvSet(dImage,cvScalar(0));
	MyPoint Orig=MLS(MyPoint(IR_X,IR_Y));
	int Orig_x=(int)(Orig.x)-(int)(pImage->width/2);
	int Orig_y=(int)(Orig.y)-(int)(pImage->height/2);

	for(i=0;i<pImage->height;i++)//遍历原图像的每个像素
	{
		for(j=0;j<pImage->width;j++)
		{
			CvScalar color;
			double x=j+IR_X;
			double y=i+IR_Y;
			MyPoint t=MLS(MyPoint(x,y));//MLS计算原图像(x,y)在目标图像的映射位置f(v)
			int m=(int)(t.x);
			int n=(int)(t.y);
			m-=Orig_x;
            n-=Orig_y;
			color=cvGet2D(pImage,i,j);//像素获取
			if(0<=m && dImage->width>m && 0<=n && dImage->height>n)
			{
			cvSet2D(dImage,n,m,color);
			}
		}
	}
图像变形算法,有正向映射和逆向映射,如果按照每个像素点,都通过上面的计算方法求取其对应变换后的像素点位置,那么其实计算量是非常大的,因为一幅图像的像素点,实在是太多了,如果每个像素点,都用上面的函数遍历过一遍,那计算量可想而知。

因此一般的变形算法是对待图像进行三角剖分:


然后只根据只对三角网格模型的顶点,根据变形算法,计算出三角网格模型每个顶点的新位置,最后再用三角形仿射变换的方法,计算三角形内每个像素点的值,得到变形后的图像,这样不仅速度快,同事解决了正向映射与逆向映射变形算法存在的不足之处,具体图像变形的正向和逆向映射存在的缺陷,可以自己查看相关的文献。

另外两种相似变换和刚性变换,可以自己查看M矩阵的计算公式,编写实现相关代码。

本文地址:http://blog.csdn.net/hjimce/article/details/46550001     作者:hjimce     联系qq:1393852684   更多资源请关注我的博客:http://blog.csdn.net/hjimce                原创文章,转载请保留本行信息*************

参考文献:

1、《Image Deformation Using Moving Least Squares》

2、3D Deformation Using Moving Least Squares

2019-05-13 20:30:13 tyfwin 阅读数 1258

 

参考资料:

陷波滤波器—matlab实现

http://blog.sina.com.cn/s/blog_ebd29d830102wdzw.html

图像复原之约束最小二乘方滤波

https://blog.csdn.net/yi_tech_blog/article/details/54605146

图像去模糊(约束最小二乘方滤波)

https://blog.csdn.net/bluecol/article/details/47359421

约束最小二乘方滤波去模糊

https://blog.csdn.net/yuyangyg/article/details/78681007

图像去模糊(维纳滤波)

https://blog.csdn.net/bluecol/article/details/46242355

 

%% 仿真逆滤波、维纳滤波,约束最小二乘滤波

close all;
clear all;
clc;
% Display the original image.
I = imread('Lena512.png'); 
[d1,d2,d3] = size(I); 
if(d3 > 1) 
I = rgb2gray(I);
end
I = im2double(I);
[hei,wid,~] = size(I);
subplot(3,3,1),imshow(I);
title('Original Image ');

% Simulate a motion blur.
LEN = 50;
THETA = 11;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');
subplot(3,3,2), imshow(blurred); title('Blurred Image');

% Simulate additive noise.
noise_mean = 0;
% noise_var = 0.00001;
noise_var = 0.0001;
blurred_noisy = imnoise(blurred, 'gaussian', ...
                        noise_mean, noise_var);
subplot(3,3,3), imshow(blurred_noisy)
title('Simulate Blur and Noise')

%% 使用自带的 deconvwnr 进行维纳滤波。 deconvreg进行约束最小二乘滤波
% deconvwnr 维纳滤波,如果没有参数NSPR,则为逆滤波

%自带 逆滤波 对已添加噪声图像 deconvwnr
deblurred4 = deconvwnr(blurred_noisy,PSF);
subplot(3,3,4), imshow(deblurred4); title('deconvwnr逆滤波 对 运动+噪声')

%自带 维纳滤波 对已添加噪声图像 deconvwnr
deblurred4 = deconvwnr(blurred_noisy,PSF,0.005); %0.005为噪声信号比
subplot(3,3,5), imshow(deblurred4); title('deconvwnr维纳滤波 对 运动+噪声')

%自带的 deconvreg 进行约束最小二乘滤波
subplot(3,3,6);
imshow(deconvreg(blurred_noisy, PSF,20)); %20 为噪声功率
title('deconvreg最小二乘滤波 对 运动+噪声');

%% 自写 逆滤波,约束最小二乘滤波

%自写 逆滤波 对未添加噪声图像
If = fft2(blurred);
Pf = psf2otf(PSF,[hei,wid]);
deblurred = ifft2(If./Pf);
subplot(3,3,7), imshow(deblurred); title('逆滤波 对 仅运动')
%自写 逆滤波 对已经添加噪声图像
If2 = fft2(blurred_noisy);
deblurred2 = ifft2(If2./Pf);
subplot(3,3,8), imshow(deblurred2); title('逆滤波 对 运动+噪声')

% Try restoration using  Home Made Constrained Least Squares Filtering.
% 自写约束最小二乘滤波
p = [0 -1 0;-1 4 -1;0 -1 0];
P = psf2otf(p,[hei,wid]);

gama = 0.001;
If = fft2(blurred_noisy);

numerator = conj(Pf);
denominator = Pf.^2 + gama*(P.^2);

deblurred2 = ifft2( numerator.*If./ denominator );
subplot(3,3,9), imshow(deblurred2)
title('约束最小二乘滤波');


 

没有更多推荐了,返回首页