卡通图像变形算法(Moving Least Squares)附源码
本文介绍一种利用移动最小二乘法来实现图像变形的方法,该方法由用户指定图像中的控制点,并通过拖拽控制点来驱动图像变形。假设p为原图像中控制点的位置,q为拖拽后控制点的位置,我们利用移动最小二乘法来为原图像上的每个像素点v构建相应的仿射变换lv(x),并通过该变换来计算得到图像变形后的位置:
其中权重wi的表达式为wi = 1/|pi - v|2α。
仿射变换lv(x)由两部分组成lv(x) = xM + T,其中M为线性转换矩阵,T为平移量。事实上将最小化表达式对变量T求偏导后可以得到T的表达式T = q* - p*M,其中p* = ∑wipi/∑wi,q* = ∑wiqi/∑wi。
于是仿射变换可以化简为lv(x) = (x - p*)M + q*,而最小化表达式可以变化为:
其中,。
注意移动最小二乘法并未对转换矩阵M进行条件限制,如果添加其他限制条件后,能得到不同形式的转换矩阵M,文章根据不同的转换矩阵M提出了三种变形方法:仿射变形(Affine Deformation)、相似变形(Similarity Deformation)和刚性变形(Rigid Deformation),下面分别介绍这三种方法。
- 仿射变形(Affine Deformation)
仿射变形是利用经典正规方程对最小化表达式直接求解得到的结果:
有了旋转矩阵M的表达式后,我们得到变形的表达式:
由于用户是通过控制q的位置来实现图像变形,而p的位置是固定不变的,因此上式中大部分内容可以预先计算并保存,从而提高运算速度,重写变形表达式如下:
其中。
% Precomputing the affine deformation: function data = Precompute_Affine(p,v,w) % Computing pstar: pstar = Precompute_pstar(p,w); % Precomputing the first matrix: M1 = v - pstar; np = size(p,1); nv = size(v,1); % Iterating on points: phat = cell(1,np); M2 = zeros(2,2,nv); for i = 1:np % Computing the hat points: phat{i} = bsxfun(@minus, p(i,:), pstar); % Computing the matrix elements: M2(1,1,:) = M2(1,1,:) + permute(w(:,i).*phat{i}(:,1).^2, [3,2,1]); M2(1,2,:) = M2(1,2,:) + permute(w(:,i).*phat{i}(:,1).*phat{i}(:,2), [3,2,1]); M2(2,1,:) = M2(1,2,:); M2(2,2,:) = M2(2,2,:) + permute(w(:,i).*phat{i}(:,2).^2, [3,2,1]); end % Computing the inverse: nv = size(v,1); IM2 = mmx('backslash', M2, repmat(eye(2), [1,1,nv])); % Computing the first product elements: F1 = [sum(M1.*squeeze(IM2(:,1,:))',2), sum(M1.*squeeze(IM2(:,2,:))',2)]; % Computing the A values: A = zeros(nv,np); for i = 1:np A(:,i) = sum(F1.*phat{i},2).*w(:,i); end % The data structure: data.A = A; end
- 相似变形(Similarity Deformation)
由于仿射变形包含非均匀缩放,因此其变形效果不是很好。相似变形是仿射变形的一个特殊子集,它的变形效果只包含平移、旋转和均匀缩放,我们限制转换矩阵M使其满足MTM = λ2I,根据该条件得到相似变形的转换矩阵M如下:
其中。
与仿射变形一样,我们提取出可以预先计算的部分后得到变形表达式:
其中。
% Precomputing the similar deformation: function data = Precompute_Similar(p,v,w) % Computing pstar: pstar = Precompute_pstar(p,w); np = size(p,1); nv = size(v,1); % Iterating on points: phat = cell(1,np); mu = zeros(nv,1); for i = 1:np % Computing the hat points: phat{i} = bsxfun(@minus, p(i,:), pstar); % Updating the values of mu: mu = mu + w(:,i).*sum(phat{i}.^2,2); end % Computing the matrix A: A = cell(1,np); R1 = v - pstar; R2 = [R1(:,2),-R1(:,1)]; for i = 1:np L1 = phat{i}; L2 = [L1(:,2),-L1(:,1)]; % [col1 col2] % [col3 col4] A{i} = [w(:,i).*sum(L1.*R1,2), ... w(:,i).*sum(L1.*R2,2), ... w(:,i).*sum(L2.*R1,2), ... w(:,i).*sum(L2.*R2,2)]; end % Premultiplying A/mu: for i = 1:np A{i} = bsxfun(@rdivide, A{i}, mu); end % The data structure: data.A = A; end
- 刚性变形(Rigid Deformation)
最近许多研究都是关于刚性变形,也就是说变形不含任何缩放效果,我们进一步限制转换矩阵M使其满足MTM = I,这样可以得到刚性变形,刚性变形的表达式如下:
其中,式中Ai表达式与相似变形中的Ai表达式相同。
% Precomputing the rigid deformation: function data = Precompute_Rigid(p,v,w) % Computing pstar: pstar = Precompute_pstar(p,w); np = size(p,1); % Iterating on points: phat = cell(1,np); for i = 1:np % Computing the hat points: phat{i} = bsxfun(@minus, p(i,:), pstar); end % Computing the matrix A: A = cell(1,np); R1 = v - pstar; R2 = [R1(:,2),-R1(:,1)]; for i = 1:np L1 = phat{i}; L2 = [L1(:,2),-L1(:,1)]; % [col1 col2] % [col3 col4] A{i} = [w(:,i).*sum(L1.*R1,2), ... w(:,i).*sum(L1.*R2,2), ... w(:,i).*sum(L2.*R1,2), ... w(:,i).*sum(L2.*R2,2)]; end % The norm of v-pstar: norm_v_pstar = sqrt(sum((v - pstar).^2, 2)); % The data structure: data.A = A; data.norm_v_pstar = norm_v_pstar; end
本文为原创,转载请注明出处:http://www.cnblogs.com/shushen
参考文献:
[1] Scott Schaefer, Travis McPhail, and Joe Warren. 2006. Image deformation using moving least squares. ACM Trans. Graph. 25, 3 (July 2006), 533-540.