本文介绍一种利用移动最小二乘法来实现图像变形的方法，该方法由用户指定图像中的控制点，并通过拖拽控制点来驱动图像变形。假设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,);
nv = size(v,);
% Iterating on points:
phat = cell(,np);
M2 = zeros(,,nv);
for i = :np
% Computing the hat points:
phat{i} = bsxfun(@minus, p(i,:), pstar); % Computing the matrix elements:
M2(,,:) = M2(,,:) + permute(w(:,i).*phat{i}(:,).^, [,,]);
M2(,,:) = M2(,,:) + permute(w(:,i).*phat{i}(:,).*phat{i}(:,), [,,]);
M2(,,:) = M2(,,:);
M2(,,:) = M2(,,:) + permute(w(:,i).*phat{i}(:,).^, [,,]);
end % Computing the inverse:
nv = size(v,);
IM2 = mmx('backslash', M2, repmat(eye(), [,,nv])); % Computing the first product elements:
F1 = [sum(M1.*squeeze(IM2(:,,:))',2), sum(M1.*squeeze(IM2(:,2,:))',)]; % Computing the A values:
A = zeros(nv,np);
for i = :np
A(:,i) = sum(F1.*phat{i},).*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,);
nv = size(v,);
% Iterating on points:
phat = cell(,np);
mu = zeros(nv,);
for i = :np
% Computing the hat points:
phat{i} = bsxfun(@minus, p(i,:), pstar); % Updating the values of mu:
mu = mu + w(:,i).*sum(phat{i}.^,);
end % Computing the matrix A:
A = cell(,np); R1 = v - pstar;
R2 = [R1(:,),-R1(:,)];
for i = :np
L1 = phat{i};
L2 = [L1(:,),-L1(:,)]; % [col1 col2]
% [col3 col4]
A{i} = [w(:,i).*sum(L1.*R1,), ...
w(:,i).*sum(L1.*R2,), ...
w(:,i).*sum(L2.*R1,), ...
w(:,i).*sum(L2.*R2,)];
end % Premultiplying A/mu:
for i = :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,);
% Iterating on points:
phat = cell(,np);
for i = :np
% Computing the hat points:
phat{i} = bsxfun(@minus, p(i,:), pstar);
end % Computing the matrix A:
A = cell(,np); R1 = v - pstar;
R2 = [R1(:,),-R1(:,)];
for i = :np
L1 = phat{i};
L2 = [L1(:,),-L1(:,)]; % [col1 col2]
% [col3 col4]
A{i} = [w(:,i).*sum(L1.*R1,), ...
w(:,i).*sum(L1.*R2,), ...
w(:,i).*sum(L2.*R1,), ...
w(:,i).*sum(L2.*R2,)];
end % The norm of v-pstar:
norm_v_pstar = sqrt(sum((v - pstar).^, )); % 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.