ARAP

gptoolbox\mesh\arap.m

先以一个三角形为例

“arap.m”

Covariance Matrix

data.CSM = covariance_scatter_matrix(ref_V,ref_F,’Energy’,energy);

In function covariance_scatter_matrix, it constructed a covarience scatter matrix :

CSM=⎛⎝⎜kx000ky000kz⎞⎠⎟T

, in which kx is a block matrix build by arap_linear_block looks like:

kx=⎛⎝⎜w12(x1−x2)+w13(x1−x3)w21(x2−x1)w31(x3−x1)w12(x1−x2)w21(x2−x1)+w23(x2−x3)w32(x3−x2)w13(x1−x3)w23(x2−x3))w31(x3−x1)+w32(x3−x2)⎞⎠⎟

For internal two dimension, the column dimension is related to the each vertex in the mesh while the row dimension is about the neighbors of each vertex. Then the covarience scatter matrix is constructed by kx,ky,kz on its diagonal which runs over x,y,z coordinates respectively. The transport is to negative the non-diagonal element. Therefore,

kxTUx=⎛⎝⎜w12(x1−x2)(x′1−x′2)+w13(x1−x3)(x′1−x′3)w21(x2−x1)(x′2−x′1)+w23(x2−x3)(x′2−x′3)w31(x3−x1)(x′3−x′1)+w32(x3−x2)(x′3−x′2)⎞⎠⎟

kxTUy=⎛⎝⎜w12(x1−x2)(y′1−y′2)+w13(x1−x3)(y′1−y′3)w21(x2−x1)(y′2−y′1)+w23(x2−x3)(y′2−y′3)w31(x3−x1)(y′3−y′1)+w32(x3−x2)(y′3−y′2)⎞⎠⎟

kxTUz=⎛⎝⎜w12(x1−x2)(z′1−z′2)+w13(x1−x3)(z′1−z′3)w21(x2−x1)(z′2−z′1)+w23(x2−x3)(z′2−z′3)w31(x3−x1)(z′3−z′1)+w32(x3−x2)(z′3−z′2)⎞⎠⎟

kyTUx=⎛⎝⎜w12(y1−y2)(x′1−x′2)+w13(y1−y3)(x′1−x′3)w21(y2−y1)(x′2−x′1)+w23(y2−y3)(x′2−x′3)w31(y3−y1)(x′3−x′1)+w32(y3−y2)(x′3−x′2)⎞⎠⎟

kyTUy=⎛⎝⎜w12(y1−y2)(y′1−y′2)+w13(y1−y3)(y′1−y′3)w21(y2−y1)(y′2−y′1)+w23(y2−y3)(y′2−y′3)w31(y3−y1)(y′3−y′1)+w32(y3−y2)(y′3−y′2)⎞⎠⎟

kyTUz=⎛⎝⎜w12(y1−y2)(z′1−z′2)+w13(y1−y3)(z′1−z′3)w21(y2−y1)(z′2−z′1)+w23(y2−y3)(z′2−z′3)w31(y3−y1)(z′3−z′1)+w32(y3−y2)(z′3−z′2)⎞⎠⎟

kzTUx=⎛⎝⎜w12(z1−z2)(x′1−x′2)+w13(z1−z3)(x′1−x′3)w21(z2−z1)(x′2−x′1)+w23(z2−z3)(x′2−x′3)w31(z3−z1)(x′3−x′1)+w32(z3−z2)(x′3−x′2)⎞⎠⎟

kzTUy=⎛⎝⎜w12(z1−z2)(y′1−y′2)+w13(z1−z3)(y′1−y′3)w21(z2−z1)(y′2−y′1)+w23(z2−z3)(y′2−y′3)w31(z3−z1)(y′3−y′1)+w32(z3−z2)(y′3−y′2)⎞⎠⎟

kzTUz=⎛⎝⎜w12(z1−z2)(z′1−z′2)+w13(z1−z3)(z′1−z′3)w21(z2−z1)(z′2−z′1)+w23(z2−z3)(z′2−z′3)w31(z3−z1)(z′3−z′1)+w32(z3−z2)(z′3−z′2)⎞⎠⎟

CSM∗U=⎛⎝⎜⎜kxTUxkyTUxkzTUxkxTUykyTUykzTUykxTUzkyTUzkzTUz⎞⎠⎟⎟

Therefore, the first dimension of CSM∗U goes through each vertex(直接下去挨着的3个). Regarding all the vertex as a group, the second dimension iterates x,y,z of eij and the third dimension traverses the x,y,z of e′ij, which are the column and row dimension of covariance matrix Si

    % compute covariance matrix elements    S = zeros(size(data.CSM,1),dim);    S(:,1:dim) = data.CSM*repmat(U,dim,1);    % dim by dim by n list of covariance matrices    SS = permute(reshape(S,[size(data.CSM,1)/dim dim dim]),[2 3 1]);

Hence, using permute function , we change the dimension going through the vertex to the last so that covariance matrix are obtained for each vertex.

Right Hand Side

Then let us look at the right hand side of formula (8)

% precompute rhs premultiplier  if ~isfield(data,'K') || isempty(data.K)    [~,data.K] = arap_rhs(ref_V,ref_F,[],'Energy',energy);    %if flat      %data.K = repdiag(ref_map,dim) * data.K;    %end  end

    %B = arap_rhs(V,F,R);    Rcol = reshape(permute(R,[3 1 2]),size(data.K,2),1);    Bcol = data.K * Rcol;    B = reshape(Bcol,[size(Bcol,1)/dim dim]);

In function arap_rhs, it constructs a matrix K,

% K #V*dim by #(F|V)*dim*dim matrix such that:
% b = K * reshape(permute(R,[3 1 2]),size(V,1)*size(V,2)*size(V,2),1);

At first, matrix K is composed of KX,KY,KZ, it is denoted as

K=⎛⎝⎜KX000KX000KXKY000KY000KYKZ000KZ000KZ⎞⎠⎟

For better explaining, we introduce an equation:

⎛⎝⎜X000X000XY000Y000YZ000Z000Z⎞⎠⎟⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜123456789⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜123456789⎞⎠⎟⎛⎝⎜XYZ⎞⎠⎟

Looking at the left hand side and substituting X with KX, the column(vertical) dimension goes through vertices in mesh, the row(horizotal) dimension travels through vertices for finding out neighbors of each vertex. For 1 in the column vector, we can extend it vertically to make the rotation elements going through each vertex.

consequently, the result is a column vector, the first dimension goes though vertices(组内), the second dimension goes through x,y,z coordinates(组间).

B = reshape(Bcol,[size(Bcol,1)/dim dim]);

This line enable the first dimension to go through vertices and the second dimension to go through x,y,z coordinates.

Minimization process

        [U,data.preF] = min_quad_with_fixed( ...          -0.5*data.L+DQ+alpha_tik*speye(size(data.L)), ...          -B+Dl,b,bc,Aeq,Beq,data.preF);

searching upwards, we find

data.L = cotmatrix(V,F);

The details of cotmatrix can be found at http://blog.csdn.net/seamanj/article/details/53774501

    DQ = sparse(size(V,1),size(V,1));    Dl = sparse(size(V,1),dim);    % Tikhonov regularization alpha    alpha_tik = 0;

Next, we solve the formula (9)

In code it uses minimizing quadratic energy method to solve the equation. Integrating both sides with respect to unknown Z, the object function is to minimize

12ZTLZ−ZTB

Note that there is a minus symbol before 0.5, that’s because the cotmatrix has negative enties on diagonal elements, which is the negative of laplacian matrix.

最后看下这个函数 arap_linear_block

这个函数充当两个角色

即可以算CSM 也可以算K,K用于线性化R, 而CSM是用来算covariance scatter matrix

K中block的转置就是CSM中的block, 为了满足CSM, 矩阵按行相加应该得到一个0行向量, 而要满足K只需要列中对应的部分相同, 则就可以让R相加

主要讲下spokes 和 spokes-and-rims 两种方法的区别
为了让block即满足CSM又要满足K, 我们将每个三角形的权值拆开, 每条边正反各迭代一次

算右边:

spokes中每条边只受两个旋转矩阵的影响, 即两个端点处旋转的影响
而spokes-and-rims 每 条边受以三个旋转矩阵的影响 .

算R:

spokes中算R1 ,对于每个三角形只需要算w12e12e′T12+w13e13e′T13,权值拆开后就是cot3e12e′T12+cot2e13e′T13

而spokes-and-rims中算R1,对于每个三角形需要算w12e12e′T12+w13e13e′T13+w23e23e′T23,权值拆开后就是cot3e12e′T12+cot2e13e′T13+cot1e23e′T23

  function K = spokes_linear_block(V,F,d)%     % Computes a matrix K such that K * R computes     %     %  鈭?wij * 0.5 * (V(i,d)-V(j,d)) * (Ri + Rj)%     % j鈭圢(i)%     % %     % Inputs:%     %   V  #V by dim list of coordinates%     %   F  #F by 3 list of triangle indices into V%     %   d  index into columns of V%     % Output:%     %   K  #V by #V matrix%     %%     E = edges(F);   %     % Build upper part of adjacency matrix where instead of a 1 for edge from i%     % to j we have the difference of position in dimension d%     A = sparse(E(:,1),E(:,2),V(E(:,1),d)-V(E(:,2),d),size(V,1),size(V,1));%     % anti-symmetric, or considers direction of edges%     A = A-A';%     % Multiply with cotangent weights (don't worry about diagonal begin wrong%     % since in A 0it's all zeros%     L = cotmatrix(V,F);%     K = L.*A;%     % correct the diagonal (notice that the sign is positive%     K = K + diag(sum(K,2));%     K = 0.5*K;%权值在cotmatrix.m已除2,这里除2是匹配8式的0.5      % triangles      C = cotangent(V,F);      i1 = F(:,1); i2 = F(:,2); i3 = F(:,3);      I = [i1;i2;i2;i3;i3;i1;i1;i2;i3];      J = [i2;i1;i3;i2;i1;i3;i1;i2;i3];      v = [ ...         C(:,3).*(V(i1,d)-V(i2,d)) ; ...         C(:,3).*(V(i2,d)-V(i1,d)); ...          C(:,1).*(V(i2,d)-V(i3,d)); ...         C(:,1).*(V(i3,d)-V(i2,d)); ...          C(:,2).*(V(i3,d)-V(i1,d)) ; ...         C(:,2).*(V(i1,d)-V(i3,d)); ...          ... % diagonal        C(:,3).*(V(i1,d)-V(i2,d)) + C(:,2).*(V(i1,d)-V(i3,d)); ...         C(:,1).*(V(i2,d)-V(i3,d)) + C(:,3).*(V(i2,d)-V(i1,d)); ...         C(:,2).*(V(i3,d)-V(i1,d)) + C(:,1).*(V(i3,d)-V(i2,d)); ...         ];      % construct and divide by 3 so laplacian can be used as is      K = sparse(I,J,v,n,n)/2;  end  function K = spokes_and_rims_linear_block(V,F,d)    % Computes a matrix K such that K * R computes    %  鈭?   -2*(cot(aij) + cot(bij) * (V(i,d)-V(j,d)) * (Ri + Rj) +    %       -2*cot(aij) * (V(i,d)-V(j,d)) * Raij +    %       -2*cot(bij) * (V(i,d)-V(j,d)) * Rbij    % j鈭圢(i)    %     % where:          vj    %              /  |  \    %             /   |   \    %            /    |    \    %           /     |     \    %          aij    |    bij    %           \     |     /    %            \    |    /    %             \fij|gij/    %              \  |  /    %                 vi    %     % Inputs:    %   V  #V by dim list of coordinates    %   F  #F by 3 list of triangle indices into V    %   d  index into columns of V    % Output:    %   K  #V by #F matrix    %    if simplex_size == 3      % triangles      C = cotangent(V,F);      i1 = F(:,1); i2 = F(:,2); i3 = F(:,3);      I = [i1;i2;i2;i3;i3;i1;i1;i2;i3];      J = [i2;i1;i3;i2;i1;i3;i1;i2;i3];      v = [ ...         C(:,3).*(V(i1,d)-V(i2,d)) + C(:,2).*(V(i1,d)-V(i3,d)); ...         -C(:,3).*(V(i1,d)-V(i2,d)) + C(:,1).*(V(i2,d)-V(i3,d)); ...          C(:,1).*(V(i2,d)-V(i3,d)) + C(:,3).*(V(i2,d)-V(i1,d)); ...         -C(:,1).*(V(i2,d)-V(i3,d)) + C(:,2).*(V(i3,d)-V(i1,d)); ...          C(:,2).*(V(i3,d)-V(i1,d)) + C(:,1).*(V(i3,d)-V(i2,d)); ...         -C(:,2).*(V(i3,d)-V(i1,d)) + C(:,3).*(V(i1,d)-V(i2,d)); ...          ... % diagonal        C(:,3).*(V(i1,d)-V(i2,d)) - C(:,2).*(V(i3,d)-V(i1,d)); ...         C(:,1).*(V(i2,d)-V(i3,d)) - C(:,3).*(V(i1,d)-V(i2,d)); ...         C(:,2).*(V(i3,d)-V(i1,d)) - C(:,1).*(V(i2,d)-V(i3,d)); ...         ];      % construct and divide by 3 so laplacian can be used as is      K = sparse(I,J,v,n,n)/3;    elseif simplex_size == 4      % tetrahedra      assert(false)    end  end