Delaunay.rar_Delaunay_Delaunay三角网_matlab坐标网_三角网_三角网 matlab

%读取数据按坐标分成两列向量 DATA = load('Dgx.dat'); X = DATA(:,3); Y = DATA(:,2); %根据规则生成边表和三角形表 x=mean(X(:)); y=mean(Y(:)); S=(Y-y).^2+(X-x).^2; S(find(S==0))=NaN; [m,n]=find(S==min(S)); Line(1,1)=m; S=(Y-Y(m,1)).^2+(X-X(m,1)).^2; S(find(S==0))=NaN; [m,n]=find(S==min(S)); Line(2,1)=m; Line(3,1)=1; SJX=[]; count=0; judge_end = []; while(1) if (ismember(0,judge_end)|| isempty(judge_end)) judge_end=[]; p=0; for ele=Line ret=0; p=p+1; ele=Line(:,p); if (ele(3)==2) continue else count=count+1; a=sqrt((X(ele(1))-X).^2+(Y(ele(1))-Y).^2); b=sqrt((X(ele(2))-X).^2+(Y(ele(2))-Y).^2); c=sqrt((X(ele(2))-X(ele(1)))^2+(Y(ele(2))-Y(ele(1)))^2); COS=(a.^2+b.^2-c^2)./(2*a.*b); [m,n]=find(COS==min(COS)); if (count>1) for sjx = SJX if (ret==1) break else if (ismember(ele(1),sjx) && ismember(ele(2),sjx)) for point=sjx' if (ret==1) break else if (point ~= ele(1) && point ~= ele(2)) k=(Y(ele(1))-Y(ele(2)))/(X(ele(1))-X(ele(2))); b=Y(ele(1))-k*X(ele(1)); flag=k*X(point)+b-Y(point); if (flag > 0) while(1) if (k*X(m)+b-Y(m) >= 0) COS(m) = NaN; if (isnan(min(COS))) ret=1; break end [m,n]=find(COS==min(COS)); else break end end end if (flag < 0) while(1) if (k*X(m)+b-Y(m) <= 0) COS(m) = NaN; if (isnan(min(COS))) ret=1; break end [m,n]=find(COS==min(COS)); else break end end end end end end end end end end if (ret == 1) judge_end = [judge_end,1]; continue else judge_end = [judge_end,0]; i=0; if (count >1) Line(3,p)=2; end count2=0; count3=0; for ele2=Line(1:2,:) i=i+1; if (isequal([ele(1);m],ele2)||isequal([m;ele(1)],ele2)) Line(3,i)=2; count2=count2+1; end if (isequal([ele(2);m],ele2)||isequal([m;ele(2)],ele2)) Line(3,i)=2; count3=count3+1; end end if (count2 == 0) Line=[Line,[ele(1);m;1]]; end if (count3 == 0) Line=[Line,[ele(2);m;1]]; end SJX=[SJX,[ele(1);ele(2);m]]; end end end else break end end %按照生成的边表画图 for pos=Line line([X(pos(1)),X(pos(2))],[Y(pos(1)),Y(pos(2))]) hold on; end