www.gusucode.com > 物流及供应链管理应用技术研发中心工具箱 > 物流及供应链管理应用技术研发中心工具箱/物流及供应链管理应用技术研发中心工具箱/LSCM2points_batch.m

    %function [x,y]=LSCM2points(v, fac, iter)


iter=1000;
%Initialise, and form face array (f). f is of form: fac(vertex(1, 2, 3),Coordinate(x, y, z), facenumber (1:length(f)));
len=size(F, 2);
verts=length(V);
f=zeros(3,3,len); 
global vec
vec=zeros(4,3,len);%Vectors as detiled below
xys=zeros(3,3,len);%xys(Vertex (1,2,3), Coordinate (x', y', z'), facenumber). Note that z should always be the same per triangle, but is retained for testing.
W=zeros(3,1);
%dT=zeros(1,len);
M=zeros(len, verts);
for face=1:len
    for point=1:3
        f(point, :, face)=V(:, F(point, face));
    end
    %Re-write of vector creation.
    vec(1,:,face)=f(3,:,face)-f(2,:,face);
    unify(1,face);
    vec(4,:,face)=f(3,:,face)-f(1,:,face);
    vec(3,:,face)=cross(vec(1,:,face), vec(4,:,face));
    unify(3,face);
    vec(2,:,face)=cross(vec(1,:,face), vec(3,:,face));
    unify(2,face);
    
    %Convert coordinate details to new orthonormal basis coordinates.
    for point=1:3
        for coor=1:3
            xys(point,coor,face)=dot(f(point,:,face),vec(coor,:,face)');%This needs correcting to fulfil theorem 7.5 of p341. Done (I think!)
        end
    end
    %Calculate W and dT as temporary numbers
    W(1)=xys(3,1,face)-xys(2,1,face)+i*(xys(3,2,face)-xys(2,2,face));
    W(2)=xys(1,1,face)-xys(3,1,face)+i*(xys(1,2,face)-xys(3,2,face));
    W(3)=xys(2,1,face)-xys(1,1,face)+i*(xys(2,2,face)-xys(1,2,face));
    dT=xys(1,1,face)*xys(2,2,face)-xys(1,2,face)*xys(2,1,face)+xys(2,1,face)*xys(3,2,face)-xys(2,2,face)*xys(3,1,face)+xys(3,1,face)*xys(1,2,face)-xys(3,2,face)*xys(1,1,face);
%    if (dT==0)
%        dT
%        W
%        face
%        xys(:,:,face)
%    end
    %Build M=mij where mij=(Wj, Ti/sqrt(dT) if vertex j belongs to triangle
    %Ti, or 0 otherwise. Rows (i) indexed by faces, cols (j) by vertices.
    for vert=1:3
        M(face,F(vert,face))=W(vert)/sqrt(dT);
    end
%    face;
end

%Choose pinned coordinates
df=diff(f(:, 3, :));
df=df(1,:,:).^2+df(2,:,:).^2;
df=squeeze(df);
pinned_face=find(df==min(df));
pinned_verts=F(:, pinned_face);
pinned_verts=pinned_verts(1:2);
%len
%size(M)
%pinned_verts=[1 verts];
pv=sort(pinned_verts);

%Decompose M into block matrices M=(Mf Mp) for free and pinned coordinates
[a,b]=size(M);
Mp=zeros(a,2,2);
Mf=zeros(a,b-2,2);
Mp(:,:,1)=real(M(:,pinned_verts));
Mp(:,:,2)=imag(M(:,pinned_verts));
Mf(:,:,1)=real(M(:,[1:pv(1)-1 pv(1)+1:pv(2)-1 pv(2)+1:end]));
Mf(:,:,2)=imag(M(:,[1:pv(1)-1 pv(1)+1:pv(2)-1 pv(2)+1:end]));

clear  W vec xys; %M not cleared for testing purposes


%Form A and b matrices (See L�vy paper for details)
Aa=cat(1,Mf(:,:,1), -Mf(:,:,2));
Ab=cat(1,Mf(:,:,2), Mf(:,:,1));
A=cat(2,Aa,Ab);
U=cat(1,f(1:2,1,round(len/2)), f(1:2,2,round(len/2)));
ba=cat(1,Mp(:,:,1), -Mp(:,:,2));
bb=cat(1,Mp(:,:,2), Mp(:,:,1));
bc=cat(2,ba,bb);
b=-bc*U;

%Calculate vector of unknowns (x1) using least squares conjugate gradient method
%flag included to suppress output message
%tic
%[x1 flag]=lsqr(A,b,[],iter);
x1 = A\b;
%toc
%x1 = pinv(A)*b;
len=length(x1);
%size(v(1,1))
%size(x1(1:len/2))
%size(v(1,end))
%x=[v(1,1) x1(1:len/2)' v(1,end)];
%y=[v(2,1) x1(1+len/2:len)' v(2,end)];
x=x1(1:len/2)';
x=[x(1:pv(1)-1) (x(pv(1)-1)+x(pv(1)))/2 x(pv(1):pv(2)-2) (x(pv(2)-2)+x(pv(2)-1))/2 x(pv(2)-1:end)];
y=x1(1+len/2:len)';
y=[y(1:pv(1)-1) (y(pv(1)-1)+y(pv(1)))/2 y(pv(1):pv(2)-2) (y(pv(2)-2)+y(pv(2)-1))/2 y(pv(2)-1:end)];

clear x1 len b bc U bb ba A Aa Ab Mp Mf M a pv pinned_verts pinned_face df vert  dT xys point coor f verts len