closecirclesearch.m

利用ＭＡＴＬＡＢ实现了对测量中的平面网

function[RT,p1,p2]=closecirclesearch(be,en,sd,gd)%优化树枝法，建立生成树，得到的RT阵
bee=be;
enn=en;
while 2
    g=size(bee,1);
    num=zeros(1,sd);
    for i=1:g%计算结点度
        num(bee(i))=num(bee(i))+1;
        num(enn(i))=num(enn(i))+1;
    end
    [degreenum,pointnum]=max(num);
    course=zeros(1,gd);%记录观测路线选择信息
    been=[bee enn];
    temp1=[];
    i1=0;
    i2=0;
    [i,j]=find(been==pointnum);
    tem=been(i,:);
    tem=reshape(tem,1,prod(size(tem)));
    tem((tem==pointnum))=[];
    temp=[pointnum tem];
    course(i)=1;
    while 1%寻找网的一棵树
        for i=tem
            [i1,j1]=find(i==been);
            beentem=been(i1,:);
            for j=1:size(beentem,1);
                xxx=beentem(j,:);
                newpoint=xxx(xxx~=i);
                if any(temp==newpoint)==0;
                    temp=[temp newpoint];
                    temp1=[temp1 newpoint];
                    course(i1(j))=1;
                end
            end
        end
        tem=temp1;
        temp1=[];
        if size(temp,2)==sd;
            break
        end
        if isempty(tem)==1;
            break
        end
    end
    if size(temp,2)==sd;
        break
    end
    [i j]=find(course==1);
    bee(i)=[];
    enn(i)=[];
end
[i j]=find(course==0);
n=length(i);
%--------------------------------------------------------------------------
G1=zeros(sd,gd-n);
G2=zeros(sd,n);
A=[];B=[];%A为树枝，B为余枝
for i=1:gd
    if course(i)==0
        B=[B;be(i) en(i)];
    else
        A=[A;be(i) en(i)];
    end
end
for i=1:gd-n%形成关联阵G1
    G1(A(i,1),i)=1;
    G1(A(i,2),i)=-1;
end
for i=1:n
    G2(B(i,1),i)=1;
    G2(B(i,2),i)=-1;
end
A=[A;B];%重新合并A,B为树枝与余枝排序后的观测信息阵A
G1(2,:)=[];
G2(2,:)=[];
RTe=-(inv(G1)*G2)';
RT=[RTe eye(size(RTe,1))];
p1=A(:,1);
p2=A(:,2);
n=size(RT);
%进行环间的运算,使得闭合环的边数最少==========
while 1
    nn=0;
    arte=abs(RT);
    nl=sum(arte,2);
    %排序
    for i=1:n-1
        for j=i+1:n
            if nl(j)>nl(i)
              xx=nl(j);
              nl(j)=nl(i);
              nl(i)=xx;
              xx=RT(j,:);
              RT(j,:)=RT(i,:);
              RT(i,:)=xx;
           end
       end
    end
    arte=abs(RT);
    for i=1:n-1
        xi=arte(i,:);
        xi((xi==0))=[];
        i0=length(xi);
        for j=i+1:n
           xj=arte(j,:);
           xj((xj==0))=[];
           j0=length(xj);
           xn=arte(i,:)-arte(j,:);
           xn((xn==0))=[];
           n0=length(xn);
           n2=(i0+j0-n0)/2;
           if n2>j0-n2
               r1=RT(i,:)+RT(j,:);
               r1((r1==0))=[];
               t1=length(r1);
               r2=RT(i,:)-RT(j,:);
               r2((r2==0))=[];
               t2=length(r2);
               if t1>t2
                   RT(i,:)=RT(i,:)-RT(j,:);
               else
                   RT(i,:)=RT(i,:)+RT(j,:);
               end
               nn=nn+1;
               break;
           end
       end
   end
   if nn==0
       break;
   end
end