erweilunkuoxian.m

对电容层析成象的ECT断层图象 进行3维显示

t = linspace(-2,2,40)*pi;
x = sin(t);
z = cos(t);
fill(x,z,'k');
axis auto
hold on
t=linspace(0,1,40)*pi;
x=sin(t);
for i=1:20;
 z(i)=-sqrt(1-x(i).^2);
end;
for i=21:40;
    z(i)=-1;
end;
fill(x,z,'w');
hold on
t=linspace(-1,0,40)*pi;
x=sin(t);
for i=1:20
z(i)=-sqrt(1-x(i).^2);
end;
for i=21:40;
    z(i)=-1;
end;
fill(x,z,'w')
hold on
t=linspace(-1,0,40)*pi
x=sin(t)
for i=1:20;
z(i)=sqrt(1-x(i).^2);
end;
for i=21:40
    z(i)=1;
end;
fill(x,z,'w');
hold on
t=linspace(0,1,40)*pi;
x=sin(t);
for i=1:20
z(i)=sqrt(1-x(i).^2);
end;
for i=21:40;
    z(i)=1;
end;
fill(x,z,'w');
axis off
hold on 
bw=bwperim(~I)
for i=1:10
    for j=1:10
           if (bw(i,j)==1)
           y=0.9-0.2*(i-1)
            x=-0.9+0.2*(j-1)
            plot(x,y,'k*')
            hold on
       end 
   end
     
end
hold on
U=zeros(10,30);
V=zeros(10,30);
Z=zeros(10,10);
F=zeros(10,10);
Q=zeros(10,10);
xt=zeros(1,2);
zt=zeros(1,2);
a=0;
b=0;
for i=1:10
    for j=1:10
            Z(i,j)=0.9-0.2*(i-1);
            X(i,j)=-0.9+0.2*(j-1);    
       end 
   end
 for i=1:10
    for j=1:10
        if(X(i,j)^2+Z(i,j)^2<=1)
            bw1(i,j)=1;
        end
    end
end 
I=~(~I&bw1);
bw=bwperim(~I);
BW=bw;        
%写坐标矩阵
k=0;
for i=1:10
    for j=1:10
         hold on 
         if(BW(i,j)==1&F(i,j)==0) 
              xb=X(i,j);
              zb=Z(i,j);
              Q(i,j)=1;
              g=0;
              quan=0;
              m=i;
              n=j;
              k=k+1;
              h=0;
              %找到起点    设起点标志位
              while quan==0
             while g==0&quan==0
             %右扫描
               xt(1)=X(i,j);
               zt(1)=Z(i,j);
               if (i~=1&i~=10&j~=10)
                   %一般情况
                  if(BW(i-1,j)==1& F(i-1,j)==0)
                    a=i-1;
                    b=j;
                  elseif(BW(i-1,j+1)==1& F(i-1,j+1)==0)
                    a=i-1;
                    b=j+1;
                  elseif(BW(i,j+1)==1& F(i,j+1)==0)
                     a=i;
                     b=j+1;
                  elseif(BW(i+1,j+1)==1&F(i+1,j+1)==0)
                     a=i+1;
                     b=j+1; 
                  elseif(BW(i+1,j)==1&F(i+1,j)==0)
                     a=i+1;
                     b=j;
                   else g=1;    
                   end          
                    %一般情况                      
               elseif i==1&j~=10
                   %特殊情况之一右扫下三点
                  if(BW(i,j+1)==1& F(i,j+1)==0)
                     a=i;
                     b=j+1;
                  elseif(BW(i+1,j+1)==1&F(i+1,j+1)==0)
                     a=i+1;
                     b=j+1; 
                  elseif(BW(i+1,j)==1&F(i+1,j)==0)
                     a=i+1;
                     b=j;
                  else   g=1;    
                   end   
                 %右扫下三点结束                 
               elseif i==10&j~=10
                   %特殊情况之二右扫上三点
                     if(BW(i-1,j)==1& F(i-1,j)==0)
                     a=i-1;
                     b=j;
                     elseif(BW(i-1,j+1)==1& F(i-1,j+1)==0)
                      a=i-1;
                     b=j+1;
                     elseif(BW(i,j+1)==1& F(i,j+1)==0)
                      a=i;
                      b=j+1;
                     else g=1;    
                     end 
                     %右扫上三点结束                     
                 elseif j==10&i~=1&i~=10
                   %特殊情况之三
                   if(BW(i-1,j)==1& F(i-1,j)==0)
                     a=i-1;
                     b=j;
                    elseif(BW(i+1,j)==1& F(i+1,j)==0)
                     a=i+1;
                     b=j;
                    else g=1;
                    %结束
                   end
               elseif j==10&i==1                   
                   if(BW(i+1,j)==1& F(i+1,j)==0)
                     a=i+1;
                     b=j;
                    else g=1;
                    end
                elseif j==10&i==10                    
                    if(BW(i-1,j)==1& F(i-1,j)==0)
                     a=i-1;
                     b=j;
                    else g=1;
                    end
                end
               if(g~=1)
                 %连线
                xt(2)=X(a,b);
                zt(2)=Z(a,b);      
                hold on
                plot(xt,zt,'r');
                hold on                
                i=a;
                j=b;
                F(i,j)=1;
                h=h+1;               
                U(k,h)=X(i,j);
                 V(k,h)=Z(i,j);         
             end
                %连线结束
                if Q(i,j)==1
                    quan=1;
                    Q(i,j)=0;
                end                
            end%  右扫描结束             
          %判别扫描是否结束                  
               if ((i~=1&i~=10&j~=1)&((BW(i+1,j-1)==1&F(i+1,j-1)==0)|(BW(i,j-1)==1&F(i,j-1)==0)|(BW(i-1,j-1)==1&F(i-1,j-1)==0)))|((i==1&j~=1)&((BW(i+1,j-1)==1&F(i+1,j-1)==0)|(BW(i,j-1)==1&F(i,j-1)==0)))|((i==10&j~=1)&((BW(i,j-1)==1&F(i,j-1)==0)|(BW(i-1,j-1)==1&F(i-1,j-1)==0))) 
                                 %左扫描  
                  while(g==1&quan==0)
                        xt(1)=X(i,j);
                        zt(1)=Z(i,j);                        
                        if (i~=1&i~=10&j~=1)
                         %一般情况
                                if(BW(i+1,j)==1&F(i+1,j)==0) 
                                   a=i+1;
                                   b=j;
                                elseif(BW(i+1,j-1)==1&F(i+1,j-1)==0)
                                   a=i+1;
                                   b=j-1;
                                elseif(BW(i,j-1)==1&F(i,j-1)==0)
                                   a=i;
                                    b=j-1;
                                elseif(BW(i-1,j-1)==1&F(i-1,j-1)==0)
                                    a=i-1;
                                    b=j-1;
                                elseif(BW(i-1,j)==1&F(i-1,j)==0)
                                    a=i-1;
                                    b=j;
                                else g=0;
                                end 
                                %一般情况                    
                        elseif i==1&j~=1
                        %特殊情况之一左扫下三点
                                 if(BW(i+1,j)==1&F(i+1,j)==0) 
                                   a=i+1;
                                   b=j;
                                 elseif(BW(i+1,j-1)==1&F(i+1,j-1)==0)
                                   a=i+1;
                                   b=j-1;
                                 elseif(BW(i,j-1)==1&F(i,j-1)==0)
                                    a=i;
                                    b=j-1;
                                  else g=0;
                                  end    
                        elseif i==10&j~=1
                       %特殊情况之二左扫上三点 
                                 if(BW(i,j-1)==1&F(i,j-1)==0)
                                   a=i;
                                    b=j-1;
                                 elseif(BW(i-1,j-1)==1&F(i-1,j-1)==0)
                                    a=i-1;
                                    b=j-1;
                                 elseif(BW(i-1,j)==1&F(i-1,j)==0)
                                    a=i-1;
                                    b=j;
                                 else g=0;
                                  end 
                        elseif j==1&i~=1&i~=10
                       %特殊情况之三
                             if(BW(i-1,j)==1& F(i-1,j)==0)
                               a=i-1;
                               b=j;
                             elseif(BW(i+1,j)==1& F(i+1,j)==0)
                               b=j;
                             else g=0;
                              end
                        elseif j==1&i==1
                            if(BW(i+1,j)==1& F(i+1,j)==0)
                              b=j;
                            else g=1;
                            end
                       elseif j==1&i==10                    
                          if(BW(i-1,j)==1& F(i-1,j)==0)
                               a=i-1;
                                  b=j;
                          else g=0;
                          end                                          
                        %结束
                    end 
                    %  左扫                                                  
               if(g~=0)
                 %连线
                xt(2)=X(a,b);
                zt(2)=Z(a,b);             
                hold on
                plot(xt,zt,'r');
                hold on               
                i=a;
                j=b;
                F(i,j)=1;
                  h=h+1;
                U(k,h)=X(i,j);
                V(k,h)=Z(i,j);                
                  if Q(i,j)==1
                    quan=1;
                    Q(i,j)=0;
                end         
                end
                %连线结束
            end                
             %判别扫描是否结束
           if ~(((i~=1&i~=10&j~=10)&((BW(i+1,j+1)==1&F(i+1,j+1)==0)|(BW(i,j+1)==1&F(i,j+1)==0)|(BW(i-1,j+1)==1&F(i-1,j+1)==0)))|((i==1&j~=10)&((BW(i+1,j+1)==1&F(i+1,j+1)==0)|(BW(i,j+1)==1&F(i,j+1)==0)))|((i==10&j~=10)&((BW(i,j+1)==1&F(i,j+1)==0)|(BW(i-1,j+1)==1&F(i-1,j+1)==0))) )
               quan=1;             
           end  
        else quan=1;               
       end                
     end %(quan==0)段程序   
    end                
  end  
end