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

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

    原文地址: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

    展开全文
  • 最小二乘法 当我们在测量某个值y时,由于误差的存在,可能多次测量的结果不尽相同 我们把多次测量得到的不同结果yi画在同一坐标系中 同时将猜测的实际值y也画在坐标系中 每个yi和y都有一个差值| y - yi |,称为误差...

    这是我的【项目笔记】利用OpenCV的MLS图像扭曲变形实现中的第一部分
    本文主要对MLS进行了一定讲解

    先简单了解一下什么是最小二乘法

    最小二乘法

    当我们在测量某个值y时,由于误差的存在,可能多次测量的结果不尽相同

    我们把多次测量得到的不同结果yi画在同一坐标系中

    同时将猜测的实际值y也画在坐标系中

    yi在坐标系中

    每个yiy都有一个差值| y - yi |称为误差

    记所有误差的平方和

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-1Pc3Ah3Q-1573289354246)(https://i.loli.net/2019/11/07/YhmLJkTn5yAIRo9.png)]

    由于实际值y是我们猜测的,所以它的值可以变化,同时误差的平方和ε也会随之改变

    y在坐标系中变化

    于是高斯或是法国科学家勒让德就提出使误差的平方和最小的 y 就是真值,这是基于,如果误差是随机的,应该围绕真值上下波动

    这就是最小二乘法,即

    最小二乘法

    此外,经证明得出误差的分布服从正态分布(不愧是天下第一分布),这里就不证明了

    总的来说,对于被选择的参数,应该使算出的函数曲线与观测值之差的平方和最小。用函数表示为:

    最小化问题的精度,依赖于所选择的函数模型

    移动最小二乘法

    移动最小二乘法与传统的最小二乘法相比,有两个比较大的改进:

    • 拟合函数的建立不同。这种方法建立拟合函数不是采用传统的多项式或其它函数,而是由一个系数向量 a(x)和基函数 p(x)构成, 这里a(x)不是常数,而是坐标 x 的函数。
    • 引入紧支( Compact Support)概念。认为点x处的值 y只受 x附近子域内节点影响,这个子域称作点 x的影响区域, 影响区域外的节点对 x的取值没有影响。在影响区域上定义一个权函数w(x),如果权函数在整个区域取为常数,就得到传统的最小二乘法。
      节选自《基于移动最小二乘法的曲线曲面拟合-曾清红》

    利用MLS变换图像

    这一部分,我主要参考了论文《Image Deformation Using Moving Least Squares》中的内容

    考虑由用户设定锚点来对图像变形进行控制的情况,首先进行准备工作,推导出公式

    准备工作

    p为一组控制点q是它对应的变形位置

    对于图像中的某一点v,有最小的仿射变换lv(x),使

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Ev2juMRP-1573289354248)(https://i.loli.net/2019/11/07/ri2nCqNtg7s9aAO.png)]

    成立。其中pi和qi是行向量,权值wi满足

    各点权重wi

    wi

    由于对于每个v都有不同的wi的值,称之为移动最小二乘最小化(a Moving Least Squares minimization)。对于每个v都有不同的lv(x)

    此时定义变形函数f(x) = lv(x),可见当v接近pi时,wi趋于无穷、f(pi) = qi,此外若qi = pi,则lv(x) = xf(v) = v

    由于lv(x)是一个仿射变换,由线性变换矩阵M和偏移量T组成,即

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-qgipx00F-1573289354250)(https://i.loli.net/2019/11/07/5TFtSdsc87HKMBi.png)]

    又可根据MT,即

    偏移量T

    T

    其中

    pq

    因此,lv(x)可以改写为

    仿射变换lv(x)

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-UTsiSVbM-1573289354255)(https://i.loli.net/2019/11/07/qwJ21uI64SY8Drt.png)]

    于是方程(1)中的最小二乘问题又可以重写为

    最小二乘

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-BHj6tzgy-1573289354259)(https://i.loli.net/2019/11/07/ISw9Ayq7Jp6b4Zo.png)]

    其中ˆpi = pi − p∗ˆqi = qi −q∗.


    可以发现,我们要进行的图像变形最终效果与变换矩阵M相关联,因此根据矩阵M的不同,可以获得不同效果的变形。论文中将其分为了仿射变换相似变换刚性变换

    在这里仅呈现了每种变换对应的变化矩阵和变化函数,具体推导可以参见论文原文

    仿射变换

    找到使方程(4)最小的矩阵M如下

    仿射变换矩阵

    M1

    由此得到的变换函数如下

    变换函数

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-simXV0Fo-1573289354261)(https://i.loli.net/2019/11/07/b6Q3OvVD42pewag.png)]

    也可以写为

    fa1_n

    其中

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-F7W3GCxo-1573289354263)(https://i.loli.net/2019/11/07/kwrxCagRHdK5qI6.png)]

    可见,一旦给定了某个点v,此处的Aj便可以提前计算得出

    如图b是仿射变换的效果图

    仿射变换效果

    相似变换

    得到进行相似变换时的变换矩阵M

    相似变换矩阵

    M2

    其中

    mius

    由此得出变换函数

    变换函数

    fa2

    其中可以提前计算的部分

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-nm6xxsRs-1573289354267)(https://i.loli.net/2019/11/09/U1TcvDLV4aYkWf7.png)]

    如图c是相似变换的效果图

    相似变换效果

    刚性变换

    刚性变换的矩阵形式与相似变化相同,仅是其中的参数miu发生了变化

    刚性变换矩阵

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-IeVTnfXA-1573289354268)(https://i.loli.net/2019/11/09/1YbZlwnScva42CB.png)]

    其中

    miur

    得出刚性变化的变化函数为

    变换函数

    fa3

    下图中d即刚性变换的效果

    刚性变换效果


    参考文献:
    《Image Deformation Using Moving Least Squares》
    《基于移动最小二乘法的曲线曲面拟合》
    参考博客:
    如何理解最小二乘法?–马同学

    展开全文
  • 在最近的项目中经常遇到给出几个点需要拟合出一条...但如果曲线曲面形式未知,可以使用移动最小二乘法或者主曲线方法。 转载:https://blog.csdn.net/liumangmao1314/article/details/54179526 MLS的讲解 移动最...

    在最近的项目中经常遇到给出几个点需要拟合出一条曲线。

    在离散的点云中,求曲线曲面拟合,不能简单地连接这些点,如果知道曲线曲面的形式,如为二次曲线等,可以简单地使用最小二乘法估计参数;但如果曲线曲面形式未知,可以使用移动最小二乘法或者主曲线方法。

    MLS的讲解

    移动最小二乘法是在最小二乘法基础上加以改进的,添加了权函数等,具体的可以参考论文,“移动最小二乘法论文”链接,这篇论文对MLS讲解的很详细,最后还给出了程序设计思路。

    我做一点点说明,论文中的矩阵A的写法欠妥,其他关于移动最小二乘法研究中还有另外一种写法:这里写图片描述,这里的B对应论文中的P,这点要注意,这样的话A就是一个矩阵。如果是线性基的二维曲线,矩阵A就是: 

    这里写图片描述,依次类推,其他的可以详看论文。

    MLS代码块

    代码的话我是根据论文中提供的程序设计,再结合一些网上的资料编写出来的,编程语言是C++;当然我也编写了python,应该是先编python,再编的C++。原因是python中可以加载一个矩阵运算库,C++中没有矩阵运算,要自己编写库,大家可以参考我这篇博客,介绍了矩阵运算链接。但是后来实验发现,python跑起来很费时间,C++只需它的一半的时间久跑完了,需要python代码也可以私信我,这里就不贴了。哦,对了,代码是包含很多自定义函数和变量,大家不要瞎贴代码,对照那篇论文的程序设计思路一下子就懂了,话不多说,上代码:

    //移动最小二乘法的具体计算过程,参照论文“基于移动最小二乘法的曲线曲面拟合”,AB矩阵参照论文“移动最小二乘法的研究”
    int MLS_Calc(int x_val,int y_val,float x[],float y[],float z[])
    {
        int max_delta=max_x-min_x;//区域半径
        float p[M][N]={0};
        float sumf[N][N]={0};
        float w[M]={0};
        for(int j=0;j<M;j++)//求w
        {
            float s=fabs((x[j]-x_val))/max_delta;
            if(s<=0.5)
                w[j]=2/3.0-4*s*s+4*s*s*s;
            else
            {
                if(s<=1)
                    w[j]=4/3.0-4*s+4*s*s-4*s*s*s/3.0;
                else
                    w[j]=0;
            }
            p[j][0]=1;//每个采样点计算基函数
            p[j][1]=x[j];
            p[j][2]=y[j];
            p[j][3]=x[j]*x[j];
            p[j][4]=x[j]*y[j];
            p[j][5]=y[j]*y[j];
        }
        f(w,x,y,sumf,p);//计算得出A矩阵
    
        float p1[N];
        Matrix A=Trans_Matrix(sumf,N);
        Matrix A_1=m_c.Matrix_copy(&A);
        m_c.Matrix_inv(&A_1);//求A矩阵的逆A_1
    
        Matrix B(N,1);//求矩阵B,N行M列
        B.init_Matrix();
        for(int j=0;j<M;j++)//求得B矩阵的每列
        {
            p1[0]=1*w[j];
            p1[1]=x[j]*w[j];
            p1[2]=y[j]*w[j];
            p1[3]=x[j]*x[j]*w[j];
            p1[4]=x[j]*y[j]*w[j];
            p1[5]=y[j]*y[j]*w[j];
            Matrix P=Trans_Matrix_One(p1,N);//数组P1转成1行N列的P矩阵
            if(j==0)//第一列直接赋值
            {
                for(int i=0;i<N;i++)
                    B.write(i,0,p1[i]);
            }
            else
            {
                m_c.Matrix_trans(&P);//矩阵转置,P转为N行1列矩阵
                m_c.Matrix_addCols(&B,&P);//矩阵B列附加,形成N行M列矩阵
            }
            P.free_Matrix();
        }
    
        float D[N]={1,x_val,y_val,x_val*x_val,x_val*y_val,y_val*y_val};
        Matrix D1=Trans_Matrix_One(D,N);//转成1行N列矩阵
    
        Matrix D_A1_mul(1,N);//定义矩阵并初始化相乘的结果矩阵,1行N列
        D_A1_mul.init_Matrix();
        if(m_c.Matrix_mul(&D1,&A_1,&D_A1_mul)==-1)
            cout<<"矩阵有误1!";//1行N列矩阵乘以N行N列矩阵得到结果为1行N列
    
        Matrix D_A1_B_mul(1,M);//定义矩阵并初始化相乘的结果矩阵,1行M列
        D_A1_B_mul.init_Matrix();
        if(m_c.Matrix_mul(&D_A1_mul,&B,&D_A1_B_mul)==-1)
            cout<<"矩阵有误2";//1行N列矩阵乘以N行M列矩阵得到记过矩阵为1行M列
    
        Matrix z1=Trans_Matrix_One(z,M);//将数组z转换成1行M列矩阵
        m_c.Matrix_trans(&z1);//转置得到M行1列矩阵
        Matrix Z(1,1);//得到矩阵结果,1行1列
        Z.init_Matrix();
        if(m_c.Matrix_mul(&D_A1_B_mul,&z1,&Z)==-1)
            cout<<"矩阵有误3!";//1行M列矩阵乘以M行1列矩阵得到1行1列矩阵,即值Z
    
        float z_val=Z.read(0,0);
        if(z_val>255)
            z_val=255;
        if(z_val<0)
            z_val=0;
    
        A.free_Matrix();
        A_1.free_Matrix();
        B.free_Matrix();
    
        D1.free_Matrix();
        D_A1_mul.free_Matrix();
        D_A1_B_mul.free_Matrix();
        z1.free_Matrix();
        Z.free_Matrix();
    
        return (int)z_val;
    }

     

     

     

     

     

    展开全文
  • 充分利用最小二乘法形函数多阶连续的特点, 能够保持在较高拟合精度的前提下获得深度图像曲面多阶连续的拟合结果。
  • 这是MLS基于C++代码实现,矩阵类运算再我的另一个资源中。代码仅供学习的,不能直接运行,需要调用。
  • 写这篇博文是应为目前为止我看到了好多领域里的经典paper算法都有涉及到移动最小二乘(MLS)。可见这个算法非常重要,先来看一下它的相关经典应用: 1、图像变形。在图像处理领域paper:《Image Deformation Using ...

    from http://www.voidcn.com/blog/hjimce/article/p-4320423.html


    一、背景意义

    写这篇博文是应为目前为止我看到了好多领域里的经典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。

    这一步实现代码如下:





    2、q*,p*的计算公式如下:

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



    3、仿射变换矩阵M的计算公式如下:

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

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



    调用方法示例:


    图像变形算法,有正向映射和逆向映射,如果按照每个像素点,都通过上面的计算方法求取其对应变换后的像素点位置,那么其实计算量是非常大的,因为一幅图像的像素点,实在是太多了,如果每个像素点,都用上面的函数遍历过一遍,那计算量可想而知。

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


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

    另外两种相似变换和刚性变换,可以自己查看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

    展开全文
  • 目录 摘要 下载代码 MLS2D 代码讲解 ... 移动最小二乘法(Moving least squares,MLS)在我另外一篇博客有介绍。这里还是主要来说说MLS怎么做曲面拟合吧。代码在mathwork中和github上都有,我所用的是ma...

    目录

     

    摘要

    下载代码

    MLS2D 代码讲解

    代码分析

    MLS2D.m

    MLS2Dshape.m

    rectangle.m

    代码改进

    原代码不足之处

    总结


    • 摘要

           移动最小二乘法(Moving least squares,MLS)在我另外一篇博客有介绍。这里还是主要来说说MLS怎么做曲面拟合吧。代码在mathwork中和github上都有,我所用的是mathwork上的代码做了部分改进。本博客在现有代码基础上做了提速和增加随机采样点拟合方式的两方面改进,使得移动最小二乘法曲面代码适用性更为广泛。

    1、提速。在i5 3450机上一幅98*144的图像以10等步长采样曲面拟合只需3s;

    2、增加随机点曲面拟合方式。代码适用一堆离散点曲面拟合。

    • 下载代码

           链接代码中分别有三种拟合计算方式:MLS1D、MLS2D和MLS3D,区别在于自变量个数不同,如自变量有一个,数据为{1,f(1);2,f(2)...},那么用MLS1D;如果自变量有两个,如{(1,1),f(1,1);(1,2),f(1,2)...},那么用MLS2D;有三个自变量,就用MLS3D。

    • MLS2D 代码讲解

    • 代码分析

           由于我的研究方向是图像处理,因此用MLS2D来做图像的曲面拟合,即给定一组整数坐标对及其对应的RGB颜色值,然后做曲面拟合。MLS2Dshape.m是具体实行代码。

    首先感谢作者辛苦编写了该代码,原代码就不贴了,对照原代码我说下我们需要修改的地方。代码主体:

    MLS2D.m

    % SET UP NODAL COORDINATES     设置采样点
    [xI,yI] = meshgrid(-2: 0.5 : 2);
    nnodes = length(xI)*length(yI);
    
    % SET UP COORDINATES OF EVALUATION POINTS  
    [x,y] =meshgrid (-2 : 0.1: 2);
    npoints = length(x)*length(y);
    scale = 3;
    % DETERMINE RADIUS OF SUPPORT OF EVERY NODE  支持半径
    dmI = scale *0.5* ones(1, nnodes);
    
    % Evaluate MLS shape function at all evaluation points x  
    [PHI, DPHIx, DPHIy] = MLS2DShape(6, nnodes, xI,yI, npoints, x,y, dmI, 'GAUSS', 3.0 ); 
    

    MLS2Dshape.m

    ......
    %%计算每个节点在设置的半径内支撑域节点权值
        DmI = dmI;
    	% DETERMINE WEIGHT FUNCTIONS AND THEIR DERIVATIVES AT EVERY NODE
    	for i = 1 : nnodes
    		[wI(i), dwdxI(i), dwdyI(i)] =rectangleWeight(type, para, 
            x(j),y(j),xI(i),yI(i),DmI(i),DmI(i));
           xII(1,i)=xI(i);
           yII(1,i)=yI(i);
        end
    ......
    %%计算基函数
    elseif (m == 6)
          p = [ones(1, nnodes); xII;yII; xII.*xII;xII.*yII;yII.*yII]; 
          pxy   = [1; x(j); y(j);x(j)*x(j);x(j)*y(j);y(j)*y(j)];
          dpdx  = [0; 1; 0;2*x(j);y(j);0];
          dpdy = [0;0;1;0;x(j);2*y(j)];
          
          B    = p .* [wI; wI; wI; wI; wI; wI];
          DBdx   = p .* [dwdxI; dwdxI;dwdxI;dwdxI; dwdxI;dwdxI];
          DBdy  = p .* [dwdyI; dwdyI;dwdyI;dwdyI; dwdyI;dwdyI];
    end
    ......
    %%构造A矩阵并求其逆
    A   = zeros (m, m);
    	DAdx  = zeros (m, m);
    	DAdy= zeros (m, m);
    	for i = 1 : nnodes
          pp = p(:,i) * p(:,i)';
          A   = A   + wI(i) * pp;
          DAdx  = DAdx  + dwdxI(i) * pp;
          DAdy = DAdy + dwdyI(i) * pp;
        end
       ARcond = rcond(A);  
    
    %% A 矩阵病态,扩大支撑域半径,并重新从前开始计算A-1
     while ARcond<=9.999999e-015  %判断条件数
           DmI=1.1*DmI;
    ......

    rectangle.m

    %%高斯核函数,很简单,忽略dwdr计算
    function [w,dwdr] = Gauss(beta,r)
    if (r>1.0)
       w     = 0.0;
       dwdr  = 0.0;
      else
       b2 = beta*beta;
       r2 = r*r;
       eb2 = exp(-b2);
    
       w     = (exp(-b2*r2) - eb2) / (1.0 - eb2);
       dwdr  = -2*b2*r*exp(-b2*r2) / (1.0 - eb2);
       
    end

    运行该代码后,可能有两个问题:

    1、rectangle函数不存在;

    2、运行后什么也没有

           作者已经单独上传了rectangle.m文件,运行后确实什么也没有,代码中也没有给显示代码。代码没有问题,我们给它稍微修改下,用来做图像的曲面拟合。代码有点长,修改了部分:

    clc
    clear all
    I=imread('22result.jpg');
    [row,col,chn]=size(I);
    % I=I(1:col,:,:);
    % 设置节点坐标
    step=10;%步长
    xII=1: step : row;
    yII=1: step : col;
    [xI,yI] = meshgrid(yII,xII);
    nnodes = size(xI,1)*size(yI,2);
    % 设置评估点的坐标
    [x,y] = meshgrid(1: 1 : col,1: 1: row);
    
    npoints = size(x,1)*size(y,2);
    scale = 30;
    % 确定每个节点的支持半径
    dmI = scale *0.5* ones(1, nnodes);
    tic
    % 评估所有评估点x的MLS形状函数
    [PHI, DPHIx, DPHIy] = MLS2DShape(6, nnodes, xI,yI, npoints, x,y, dmI, 'GAUSS', 3.0 ); 
    toc
     
    % 曲线拟合. y = peaks(x,y)
    ZII  =I(xII,yII,:);    % 节点函数值 
    % z  =x.*exp(-x.^2- y.^2);% 确切的解决方案
    Zpoints=zeros(1,npoints);
    xh=zeros(1,npoints);
    yh=zeros(1,npoints);
    II=I-I;
    for i=1:npoints
    % Zpoints(1,i)=z(i);
    xh(1,i)=x(i);
    yh(1,i)=y(i);
    end
    Znodes=zeros(1,nnodes);
    for j=1:1
        ZI=ZII(:,:,j);
        for i=1:nnodes
            Znodes(1,i)=ZI(i);
        end                        %将二维数据转换为一维数据
        zh = PHI *Znodes';  % 逼近函数
        II(:,:,j)=reshape(zh,row,col);
    end
    ZI=double(ZI);
    plot3( xI, yI, ZI,'k.','LineWidth',2);
    hold on
    surf(x,y,II(:,:,j));
    toc

    MLS曲面拟合结果图,10等步长采样:

    运行后发现一个小问题:运行时间很长。为此,我分析作者写的代码有些耗时之处:

    1、MLS2Dshape.m中计算每个节点的支撑域权值w是用循环;

    2、MLS2Dshape.m中当矩阵A是病态时,程序会重新计算基函数;

    3、MLS2Dshape.m中用循环方法计算A。

    我也就在此基础上进行了修改,可能还有其他优化的地方。

    • 代码改进

    1、计算支撑域权值w用矩阵计算;

    2、基函数只计算一次;

    3、矩阵A的计算替换循环。

    4、当矩阵A病态时,增大它的扩大半径系数。

    详细代码我上传到我的资源里了,大家可以去下载运行,在i5 3450机上一个98*144的图片以10个等步采样点拟合只需要3s。

    • 原代码不足之处

           从原代码可以看到,代码作者只是做了等步长采样作为节点,但是有时候我们得到了一堆乱七八糟的散列点,或者图像曲面拟合时并不是方阵规则的区域,也想做曲面拟合,这时就需要对代码进行改进。限于篇幅,这里就不放了,在我资源里有,大家可以去下载,如果没有C币或积分,也可以发邮件直接问我要。不规则区域曲面拟合和随机采样点拟合类似,在不规则区域随机采样后即可拟合。代码下载链接:https://download.csdn.net/download/liumangmao1314/11132986


    2019-05-21 更新

    代码上传了百度网盘,有积分的还是给点积分吧。   仅供学习参考,禁止商业交易,请尊重作者~

    提速代码 链接: https://pan.baidu.com/s/1GUmEK0AsZdnDcOT5ZaO50Q 提取码: b3d4 

    随机采点 链接: https://pan.baidu.com/s/1G2c0keqs32RDMt72-CkcaA       提取码: k4ii 


    运行结果图如下(200个随机采样点):

    • 总结

           本博客首先分析了移动最小二乘法曲面拟合的源代码MLS2D,并以源代码MLS2D的基础上进行了提速和拟合方式改进。改进后,速度上有较大的提升;同时,拟合方式增加一种随机采样点方式使得代码满足离散点拟合和不规则区域曲面拟合。

    ---------------------

    热忱欢迎对图像矢量化有兴趣的同学一起沟通交流。

    展开全文
  • 引言陆陆续续在计算摄影学接触了不少保边滤波器,其重要性自不必说,可以用在图像的增强,图像抽象画,高动态范围图像压缩,图像色调映射等。
  • 这里要跟大家分享的paper为基于特征线的图像 morphing,对应的英文文献为《Feature-Based Image Metamorphosis》,是1992年SIGGRAPH 上的一篇paper,比较老的一篇paper,然而这篇paper引用率非常高,用于图像变形...
  • 什么是图像噪声2. 图像噪声来源3. 噪声分类及其消除方法3.1 高斯噪声及其消除方式3.1.1 如何确定一种噪声是高斯噪声3.1.2 高斯滤波消除高斯噪声3.1.3 高斯函数的重要性质3.2 瑞利噪声及其消除方式3.3 伽马(爱尔兰)...
  • 这一部分主要是最后的Problem比较重要。 带条件的代价函数ConditionedCostFunction 这个类使用户可以在使用封装好的代价函数的同时,对残差值加入一定的条件限制。举个例子,现在已经有一个代价函数可以产生N个值...
  • 本文在撰写过程中参考了由何东健教授主编、西安电子科技大学出版社出版的《数字图像处理》(第三版),一切著作权归原书作者和出版社所有。特别感谢长安大学软件系老师的认真负责的教导。 第1章 概论 1.1 数字...
  • 原文...传统图像处理部分 图像处理基础知识 彩色图像、灰度图像、二...
  • 基于Matlab的虹膜识别系统的关键技术(2) 虹膜图像的定位主要包括虹膜内边界和外边界的定位。人眼图像中瞳孔的灰度值最小,巩膜的灰度值最大,虹膜的灰度值介于二者之间。...实验中,采取基于最小二乘法原...
  • 由单目传感器获取场景图像,利用颜色信息提取路径,采用最小二乘法拟合路径参数,简化图像处理过程,提高了算法的实时性。通过消除相对参考路径的距离偏差和角度偏差来修正机器人的位姿状态,实现机器人对路径的跟踪...
  • 1、数字图像: 数字图像,又称为数码图像或数位图像,是二维图像用有限数字数值像素的表示。数字图像是由模拟图像数字化...2、数字图像处理包括内容: 图像数字化;图像变换;图像增强;图像恢复;图像压缩编码...
  • 图像细节增强

    2017-09-22 21:56:47
    主要有:L0平滑,引导滤波,快速双边滤波和边缘保留多尺度图像分解的新方法(基于加权最小二乘法框架)。 图像细节增强的核心是将原图像表示为基本分量(base layer)与细节分量(detail layer)之和,在此基础上...
  • 数字图像处理第四章数字图像处理---图像复原(三)仅有噪声的复原——空间滤波3.1 空间噪声滤波器4.3.2 自适应空间滤波 数字图像处理—图像复原 (三)仅有噪声的复原——空间滤波 如果出现的退化仅仅是噪声,...
  • 1、数字图像: ...2、数字图像处理包括内容: 图像数字化;图像变换;图像增强;图像恢复;图像压缩编码;图像分割;图像分析与描述;图像的识别分类。 3、数字图像处理系统包括部分: 输入(采集);存储
1 2 3 4 5 ... 20
收藏数 2,012
精华内容 804