www.gusucode.com > 物流及供应链管理应用技术研发中心工具箱 > 物流及供应链管理应用技术研发中心工具箱/物流及供应链管理应用技术研发中心工具箱/LSCM2points.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(:, fac(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=fac(:, 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