图像处理（十）基于特征线的图像变形-Siggraph 1992

这里要跟大家分享的paper为基于特征线的图像 morphing，对应的英文文献为《Feature-Based Image Metamorphosis》，是1992年SIGGRAPH 上的一篇paper，比较老的一篇paper，然而这篇paper引用率非常高，用于图像变形效果还是挺不错的，这个算法一般用于图像的morphing。因为这篇paper算法原理简单，易于实现，所以不用怕学习这个算法需要多长的时间。

开始之前先声明一下，这篇博文主要参考自：http://www.csie.ntu.edu.tw/~b97074/vfx_html/hw1.html  同时结合我自己的理解，跟大家分享算法，帮助更多爱好者学习，非商业用途。

一、相关理论

我们知道图像变形的本质，其实就是求取光流场，就是求取目标图像的每一个像素点在原图像的对应位置点，然后同过双线性插值的方法就可以求得目标图像。

1、单线段约束变形：

给定原图像Source Image，我们希望把原图像上的像素点P'、Q'位置移动到P、Q位置，那么其它的像素点的位置要怎么移动，才能使得得到的结果图像Destination Image不会发生严重扭曲，这便是图像变形的研究内容。如图所示，如果用逆向映射的变形方法，对于目标图像上的任意一点X，我们只需要求取source image 上的对应点X'就可以了，当然逆向映射X'往往不是整数，需要经过线性插值，获取X'的像素值。

已知PQ，P'Q'，X点，我们要怎么求出X的对应点X'呢？

其实很简单，其原理是通过保证图中二维(u,v)坐标不变就可以了，也就是我们可以通过上面的三个计算公式，求出X'。说的简单一点呢，就是保证X相对于PQ的比例位置坐标(u,v)不变。

2、多线段约束变形

上面是对于单线段约束而言的，对于多线段的情况，主要是通过加权平均的方法。

如图，现在已知P1Q1，P2Q2，X，以及源图像的P1'Q1',P2'Q2'，我们要求取X'点。

这个时候我们可以先用单线段约束的方法

(1)通过P1Q1 、P1'Q1'、X 计算出X1’;

(2)通过P2Q2 、P2'Q2'、X 计算出X2’;

然后通过加权平均的方法，求出X'：

其中权值w的计算方法就是通过点X到线段的距离成反比的函数计算：

其中length表示线段的长度，dist表示点X到线段的最短距离。a，b，p为常数，对于它们的取值我们可以选p = 0 , a = 1 , b = 2。

ok，到了这里算法就结束了，感觉松松，我们就可以计算出X’点，

因为X'计算出来一般不可能刚好位于原图像像素点的位置，因此我们需要通过双线性插值的方法，求取X'的像素值。

看完上面应该知道怎么计算X'点了吧。接着我们要进入算法实现阶段。

二、算法实现

说明以下代码参考自：http://www.csie.ntu.edu.tw/~b97074/vfx_html/hw1.html

Algorithm:

1、根据公式1,2，计算X点相对于各线段的位置(u,v)坐标。

[cpp] view plain copy

1. double new_u = dst_line.Getu(X);  
2. double new_v = dst_line.Getv(X);  

[cpp] view plain copy

1. //公式1  
2. double Line::Getu(Vector2 X)  
3. {  
4.     double X_P_x = X.x - P.x;   
5.     double X_P_y = X.y - P.y;  
6.     double Q_P_x = Q.x - P.x;  
7.     double Q_P_y = Q.y - P.y;  
8.     double u = ((X_P_x * Q_P_x) + (X_P_y * Q_P_y)) / (len * len)  ;  
9.     return u ;  
10. }  
11. //公式2  
12. double Line::Getv(Vector2 X){  
13.     double X_P_x = X.x - P.x;  
14.     double X_P_y = X.y - P.y;  
15.     double Q_P_x = Q.x - P.x;  
16.     double Q_P_y = Q.y - P.y;  
17.     double Perp_Q_P_x = Q_P_y ;    
18.     double Perp_Q_P_y = -Q_P_x ;  
19.     double v = ((X_P_x * Perp_Q_P_x) + (X_P_y * Perp_Q_P_y))/len ;   
20.     return v ;   
21. }  

2、根据公式3，然后反算X点在源图像的位置对应点X'。

[cpp] view plain copy

1. Vector2 src_point = src_line.Get_Point(new_u , new_v);  

[cpp] view plain copy

1. //根据u,v坐标可计算出pq线段的对应点x  
2. Vector2 Line::Get_Point(double u , double v)  
3. {  
4.     double Q_P_x = Q.x - P.x;  
5.     double Q_P_y = Q.y - P.y;  
6.     double Perp_Q_P_x = Q_P_y ;    
7.     double Perp_Q_P_y = -Q_P_x ;  
8.     double Point_x = P.x + u * (Q.x - P.x) + ((v * Perp_Q_P_x)/len) ;  
9.     double Point_y = P.y + u * (Q.y - P.y) + ((v * Perp_Q_P_y)/len) ;  
10.     Vector2 X;  
11.     X.x = Point_x;  
12.     X.y = Point_y;  
13.     return X ;  
14. }  

3、计算各个线段的权重。

[cpp] view plain copy

1. double src_weight = dst_line.Get_Weight(dst_point);  

[cpp] view plain copy

1. double Line::Get_Weight(Vector2 X )  
2. {  
3.     double a = parameter_a;  
4.     double b = parameter_b;  
5.     double p = parameter_p;  
6.     double d = 0.0;  
7.   
8.     double u = Getu(X);  
9.     if(u > 1.0 )  
10.         d = sqrt((X.x - Q.x) * (X.x - Q.x) + (X.y - Q.y) * (X.y - Q.y));  
11.     else if(u < 0)  
12.         d = sqrt((X.x - P.x) * (X.x - P.x) + (X.y - P.y) * (X.y - P.y));  
13.     else  
14.         d = abs(Getv(X));  
15.   
16.   
17.     double weight =pow(pow((float)len,(float)p)/(a + d) , b);  
18.     return weight;   
19. }  

然后对所有的X'点进行加权求和就可以了。

ok，上面过程的代码合在一起，遍历每一条约束线段。

4、最后进行双线性插值。

双线性插值函数如下：

[cpp] view plain copy

1. void bilinear(BitmapData *psrcImgData,float X ,float Y,byte*resultpiexl)  
2. {  
3.     int x_floor = (int)X ;  
4.     int y_floor = (int)Y ;  
5.     int x_ceil = x_floor + 1 ;  
6.     int y_ceil = y_floor + 1 ;  
7.     float a = X - x_floor ;  
8.     float b = Y - y_floor ;  
9.   
10.     if(x_ceil >= psrcImgData->Width-1)   
11.         x_ceil =psrcImgData->Width-1 ;  
12.     if(y_ceil >= psrcImgData->Height-1)   
13.         y_ceil = psrcImgData->Height-1 ;  
14.     byte leftdown[3];   
15.     byte lefttop[3];   
16.     byte rightdown[3];   
17.     byte righttop[3];   
18.     Get2D(psrcImgData,y_floor,x_floor,leftdown);  
19.     Get2D(psrcImgData,y_ceil,x_floor,lefttop);  
20.     Get2D(psrcImgData,y_floor,x_ceil,rightdown);  
21.     Get2D(psrcImgData,y_ceil,x_ceil,righttop);  
22.     for(int i = 0 ; i < 3 ; i ++)  
23.     {  
24.         resultpiexl[i] = (1-a)*(1-b)*leftdown[i] + a*(1-b)*rightdown[i] + a*b*righttop[i] + (1-a)*b*lefttop[i];  
25.     }  
26. }  
27. void Get2D(BitmapData *psrcImgData, int Y,int X, byte*piexl)  
28. {  
29.     byte*pdata=(byte*)psrcImgData->Scan0+(psrcImgData->Width*Y+X)*4;  
30.     for (int i=0;i<3;i++)  
31.     {  
32.         piexl[i]=pdata[i];  
33.     }  
34.   
35. }  

最后贴一下整个过程的代码：

[cpp] view plain copy

1. int nWidth=prightImgData->Width;  
2. int nHeight=prightImgData->Height;  
3.    for(int x = 0 ; x < nWidth ; x++)  
4. {  
5.        for(int y = 0 ; y < nHeight ; y++)  
6.     {  
7.         Vector2 dst_point ;  
8.         dst_point.x= x ;   
9.         dst_point.y= y;  
10.         double leftXSum_x = 0.0;  
11.         double leftXSum_y = 0.0;  
12.         double leftWeightSum = 0.0;  
13.         double rightXSum_x = 0.0;  
14.         double rightXSum_y = 0.0;  
15.         double rightWeightSum = 0.0;  
16.         for(int i = 0 ; i < pairs.size() ; i++)  
17.         {  
18.   
19.             Line src_line = pairs[i].leftLine;    
20.             Line dst_line = pairs[i].rightLine;  
21.   
22.             double new_u = dst_line.Getu(X);//计算（u，v）坐标  
23.             double new_v = dst_line.Getv(X);  
24.   
25.             Vector2 src_point = src_line.Get_Point(new_u , new_v);//计算源图像的对应点X'  
26.             double src_weight = dst_line.Get_Weight(dst_point);//计算权重  
27.             leftXSum_x = leftXSum_x + (double)src_point.x * src_weight ;//加权求X'的平均位置  
28.             leftXSum_y = leftXSum_y + (double)src_point.y * src_weight ;  
29.             leftWeightSum = leftWeightSum + src_weight ;  
30.         }  
31.         double left_src_x = leftXSum_x / leftWeightSum;  
32.         double left_src_y = leftXSum_y / leftWeightSum;  
33.         double right_src_x = x;  
34.         double right_src_y = y;  
35.   
36.   
37.         if(left_src_x<0)//判断是否越界  
38.             left_src_x=0;  
39.         if(left_src_y<0)  
40.             left_src_y=0;  
41.         if(left_src_x>=pleftImgData->Width)  
42.             left_src_x=pleftImgData->Width-1;  
43.         if(left_src_y>=pleftImgData->Height)  
44.             left_src_y=pleftImgData->Height-1;  
45.   
46.         byte leftimg[3];//存储最后的(x,y)点的像素值  
47.         bilinear(pleftImgData,left_src_x,left_src_y,leftimg);//线性插值  
48.         for (int i=0;i<3;i++)  
49.         {  
50.             float newpiexl=leftimg[i];  
51.   
52.   
53.         }  
54.        }  
55.    }  

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

最后看一下，用这个算法实现的变形融合：

原图像：

变形融合结果：

参考文献：

1、http://www.csie.ntu.edu.tw/~b97074/vfx_html/hw1.html
2、《Feature-Based Image Metamorphosis》