隨著社會的快速發展，我們的生活越來越離不開手機，手機已成為我們生活的必需品。大家有沒有遇到這樣的問題？就是在高樓大廈林立的地方，我們會發現手機信號特別差，通話斷斷續續并伴有雜音，視頻卡成狗，地下車庫更是盲區，為什么會出現這樣的問題呢？

其實手機信號的傳輸是一種無線電的發射與接收過程，傳播路徑主要是以直線傳播為主，若在傳播過程中遇到阻礙（即與高樓大廈墻壁發生多次反射、繞射等物理現象。其中，繞射是指在電磁波傳播路徑上，當電波被尺寸較大(與波長相比)的障礙物遮擋時，電磁波改變傳播方向的現象），就會大大削弱其信號強度。那么，該如何解決呢？

在現實生活中，運營商主要通過安裝基站的方式增強手機信號。那么問題來了，如何最優選擇基站的位置將成為重點考慮的問題？

目前，比較有代表性的就是射線跟蹤模型。射線跟蹤是一種被廣泛用于移動通信中的預測無線電波傳播特性的技術，由于移動通信中使用的超高頻微波和光同屬電磁波，有一定近似性，按光學方法辨認出多路徑信道中收、發射機間所有主要的傳播路徑。一旦這些傳播路徑被辨認后，就可根據電波傳播理論來計算每條傳播路徑信號的幅度、相位、延遲和極化，然后結合天線方向圖和系統帶寬就可得到到達接收點的所有傳播路徑的相干合成結果。

對于城市基站的二維模型，建筑群可被劃分為一定的“塊”，建筑物(即下圖中帶有灰色陰影的多邊形)則被定義為“多邊形”，多邊形的“邊”代表建筑物的表面，多邊形的“頂點”則代表了建筑物的拐角。這種簡化了的市區平面圖大致反映出城市的主體結構，利用它進行射線跟蹤，可以得到較為準確的路徑損耗。

在多邊形的頂點上僅能產生繞射，而在多邊形的邊上僅能產生反射，這些多次的反射、繞射及其組合便是收、發射機間的傳播路徑。二維射線跟蹤模型可以通過以下兩種規律分別確定反射傳播路徑和繞射傳播路徑：

(1)反射傳播路徑，如下圖(a)所示，產生反射時入射角 等于反射角 ;

(2)繞射傳播路徑，如下圖(b)所示，不論入射線以任意角度入射到建筑物頂點上，繞射射線都會以任意出射角向沒有建筑物覆蓋的區域傳播。

現在我們要求出，在發射機Tx 即坐標為(500, 200)、接收機Rx 即坐標為(250, 350)之間，通過多次反射與繞射的組合，從而形成的主要傳播路徑。

拿到這樣的題目，我們通過“枚舉法”，即發射機Tx，從0°逐步加上步長（比如0.1°，甚至可以更小，步長越小精度越高，但是程序運行時間更長），從這個角度發出一條射線，然后考慮是反射還是繞射情況，直到最后的射線與接收機Rx非常接近（可以利用點到直線的距離公式求得），此次循環才能結束，進入下一個循環，最終當Tx角度超過360°，程序才能結束。

下面，我們通過Matlab編程，運行主程序main.m文件，由于該程序涉及到大量矩陣求逆運算等過程，計算量大，在筆者電腦上面運行該程序大約需要15分鐘（取決于電腦配置），執行過程中請大家耐心等待哈。比如，點擊運行程序后，可以先去洗個澡，回來后就能看到運行結果啦

最后，我們得到運行結果為：

詳細代碼如下：

（1）主程序main.m詳細代碼

clc;
clear;
close all;


figure
building={ [5.54892205638474e+002    1.53797468354430e+002;
            4.03648424543947e+002    1.53797468354430e+002;
            4.03648424543947e+002    1.00632911392405e+002;
            5.54892205638474e+002    1.00632911392405e+002;
            5.54892205638474e+002    1.53797468354430e+002];

            [5.54892205638474e+002    1.75738396624473e+002;
            5.54892205638474e+002    1.90928270042194e+002;
            5.29684908789386e+002    1.90928270042194e+002;
            5.29684908789386e+002    1.75738396624473e+002;
            5.54892205638474e+002    1.75738396624473e+002];

            [5.54892205638474e+002    1.95147679324895e+002;
            5.54892205638474e+002    2.59282700421941e+002;
            5.36981757877280e+002    2.59282700421941e+002;
            5.36981757877280e+002    1.95147679324895e+002;
            5.54892205638474e+002    1.95147679324895e+002];

            [4.42786069651741e+002    1.75738396624473e+002;
            4.42786069651741e+002    2.59282700421941e+002;
            4.03648424543947e+002    2.59282700421941e+002;
            4.03648424543947e+002    1.75738396624473e+002;
            4.42786069651741e+002    1.75738396624473e+002];

            [5.54892205638474e+002    2.81223628691983e+002;
            5.54892205638474e+002    3.00632911392405e+002;
            4.99170812603648e+002    3.00632911392405e+002;      %%%%%%淇敼
            4.99170812603648e+002    2.81223628691983e+002;
            5.54892205638474e+002    2.81223628691983e+002];

            [4.89220563847430e+002    2.81223628691983e+002;
            4.89220563847430e+002    2.95569620253165e+002;
            4.56053067993366e+002    2.95569620253165e+002;
            4.56053067993366e+002    2.81223628691983e+002;
            4.89220563847430e+002    2.81223628691983e+002];

            [4.32172470978441e+002    2.81223628691983e+002;
            4.32172470978441e+002    3.39451476793249e+002;
            4.03648424543947e+002    3.39451476793249e+002;
            4.03648424543947e+002    2.81223628691983e+002;
            4.32172470978441e+002    2.81223628691983e+002];

            [5.54892205638474e+002    3.09071729957806e+002;
            5.54892205638474e+002    3.39451476793249e+002;
            4.47429519071310e+002    3.39451476793249e+002;
            4.47429519071310e+002    3.09071729957806e+002;
            5.54892205638474e+002    3.09071729957806e+002];

            [4.75953565505804e+002    3.62236286919831e+002;
            4.75953565505804e+002    3.79113924050633e+002;
            4.71310116086236e+002    3.79113924050633e+002;
            4.71310116086236e+002    3.62236286919831e+002;
            4.75953565505804e+002    3.62236286919831e+002];

            [4.62686567164179e+002    3.62236286919831e+002;
            4.62686567164179e+002    4.18776371308017e+002;
            4.01658374792703e+002    4.18776371308017e+002;
            4.01658374792703e+002    3.62236286919831e+002;
            4.62686567164179e+002    3.62236286919831e+002];

            [4.99834162520730e+002    3.88396624472574e+002;
            4.99834162520730e+002    4.18776371308017e+002;
            4.80597014925373e+002    4.18776371308017e+002;
            4.80597014925373e+002    3.88396624472574e+002;
            4.99834162520730e+002    3.88396624472574e+002];

            [3.79104477611940e+002    4.18776371308017e+002;
            2.94195688225539e+002    4.18776371308017e+002;
            2.94195688225539e+002    3.62236286919831e+002;
            3.79104477611940e+002    3.62236286919831e+002;
            3.79104477611940e+002    4.18776371308017e+002];

            [3.79104477611940e+002    2.81223628691983e+002;
            3.79104477611940e+002    3.39451476793249e+002;
            2.12603648424544e+002    3.39451476793249e+002;
            2.12603648424544e+002    2.81223628691983e+002;
            3.79104477611940e+002    2.81223628691983e+002];

            [3.79104477611940e+002    1.53797468354430e+002;
            3.51243781094527e+002    1.53797468354430e+002;
            3.51243781094527e+002    1.00632911392405e+002;
            3.79104477611940e+002    1.00632911392405e+002;
            3.79104477611940e+002    1.53797468354430e+002];

            [3.79104477611940e+002    1.75738396624473e+002;
            3.79104477611940e+002    2.59282700421941e+002;
            3.32669983416252e+002    2.59282700421941e+002;
            3.32669983416252e+002    2.16244725738397e+002;
            3.57213930348259e+002    2.16244725738397e+002;
            3.57213930348259e+002    1.75738396624473e+002;
            3.79104477611940e+002    1.75738396624473e+002];

            [3.38640132669983e+002    1.53797468354430e+002;
            2.41791044776119e+002    1.53797468354430e+002;
            2.41791044776119e+002    1.25949367088608e+002;
            3.02819237147595e+002    1.25949367088608e+002;
            3.02819237147595e+002    1.00632911392405e+002;
            3.27363184079602e+002    1.00632911392405e+002;
            3.27363184079602e+002    1.25949367088608e+002;
            3.38640132669983e+002    1.25949367088608e+002;
            3.38640132669983e+002    1.53797468354430e+002];

            [2.90215588723051e+002    1.75738396624473e+002;
            2.90215588723051e+002    2.08649789029536e+002;
            2.65008291873964e+002    2.08649789029536e+002;
            2.65008291873964e+002    2.59282700421941e+002;
            2.12603648424544e+002    2.59282700421941e+002;
            2.12603648424544e+002    2.33966244725738e+002;
            2.43117744610282e+002    2.33966244725738e+002;
            2.43117744610282e+002    2.01054852320675e+002;
            2.12603648424544e+002    2.01054852320675e+002;
            2.12603648424544e+002    1.75738396624473e+002;
            2.90215588723051e+002    1.75738396624473e+002];

            };
N=17;


for i=1:N
    plot(building{i}(:,1),building{i}(:,2));
    fill(building{i}(:,1),building{i}(:,2),[0.5 0.5 0.5],'FaceAlpha',0.8);
    hold on
end
axis equal
xlabel('X(m)');
ylabel('Y(m)');
%grid on


TX=[500;200];
plot(TX(1),TX(2),'*r');
h_text1 = text(TX(1)-6,TX(2)-10,['Tx'],'color','r');
set(h_text1,'fontsize',14)    %設置字體大小


RX=[250;350];
plot(RX(1),RX(2),'or');     
h_text2 = text(RX(1)-5,RX(2)+10,['Rx'],'color','r');
set(h_text2,'fontsize',14)    %設置字體大小    


[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num
    [bb]=calculate_point(aa(i,1),aa(i,2),building);
    num2=size(bb);
    number=num2(1);
    for j=1:number
         for thea=0:0.1:2*pi
            x1=bb(j,1);y1=bb(j,2);x2=x1+cos(thea);y2=y1+sin(thea); 
            node=calculate_node(x1,y1,building);
            if x2<=max(building{node}(:,1))&&x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
               [insect1,insect2]=insection_point(x1,y1,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   continue
               end
               value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
               zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
               if zhi >=0.95&&zhi<=1&&value1==1
                   plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'k');
                   plot([aa(i,1),bb(j,1)],[aa(i,2),bb(j,2)],'k');
                   plot([bb(j,1),insect1(1)],[bb(j,2),insect1(2)],'k');
                   plot([insect1(1),RX(1)],[insect1(2),RX(2)],'k');
               end
        end
    end     
end


for thea=0:0.1:2*pi
    x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
    [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
    if insect1(1)==0&&insect1(2)==0
        continue
    end
    x1=insect1(1);y1=insect1(2);
    x2=insect2(1);y2=insect2(2);    
    [x,y]=line_to_point(x1,y1,x2,y2,building);
    if x==0&&y==0
        continue
    end
    [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
    if x1==0&&y1==0
        continue
    end
   [bb]=calculate_point(x,y,building);
    num2=size(bb);
    number=num2(1);
    for j=1:number     
        value1=calculate(bb(j,1),bb(j,2),RX(1),RX(2),building);
        if value1==1
            plot([TX(1),x1],[TX(2),y1],'r');
            plot([x1,x],[y1,y],'r');
            plot([x,bb(j,1)],[y,bb(j,2)],'r');
            plot([bb(j,1),RX(1)],[bb(j,2),RX(2)],'r');
        end
    end  
end        


[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num
    for thea=0:0.1:2*pi
        x1=aa(i,1);y1=aa(i,2);x2=x1+cos(thea);y2=y1+sin(thea);     
        node=calculate_node(x1,y1,building);
        if x2<=max(building{node}(:,1))&&x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
            continue
        end
        [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
        if insect1(1)==0&&insect1(2)==0
            continue
        end
        x1=insect1(1);y1=insect1(2);
        x2=insect2(1);y2=insect2(2);    
        [x,y]=line_to_point(x1,y1,x2,y2,building);
        if x==0&&y==0
            continue
        end
        [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
        if x1==0&&y1==0
            continue
        end  
        value1=calculate(x,y,RX(1),RX(2),building);
        if value1==1
            plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'r');
            plot([aa(i,1),x1],[aa(i,2),y1],'r');
            plot([x1,x],[y1,y],'r');
            plot([x,RX(1)],[y,RX(2)],'r');
        end
    end  
end

[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num
    [bb]=calculate_point(aa(i,1),aa(i,2),building); 
    num2=size(bb);
    number=num2(1);
    for j=1:number
        value1=calculate(bb(j,1),bb(j,2),RX(1),RX(2),building);
        if value1==1
            plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'g');
            plot([aa(i,1),bb(j,1)],[aa(i,2),bb(j,2)],'g');
            plot([bb(j,1),RX(1)],[bb(j,2),RX(2)],'g');
        end
    end
end      


[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num  
    for thea=0:0.1:2*pi
         for k=1:4
            x1=aa(i,1);y1=aa(i,2);
            x2=x1+cos(thea);y2=y1+sin(thea);   
            node=calculate_node(x1,y1,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x1,y1,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x1=insect1(1);y1=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'k');
               xx=aa(i,1);yy=aa(i,2);
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end
end


 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
         for k=1:3

             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);    
            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa1=x1;bb1=y1;
            aa2=x;bb2=y;
            x2=x+cos(thea1);y2=y+sin(thea1);   
            node=calculate_node(x,y,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x,y,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x=insect1(1);y=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa1],[TX(2),bb1],'k');
               plot([aa1,aa2],[bb1,bb2],'k');
               xx=aa2;yy=bb2;
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end    
 end           


 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
         for k=1:2

             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa3=x1;bb3=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 

            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa1=x1;bb1=y1;
            aa2=x;bb2=y;
            x2=x+cos(thea1);y2=y+sin(thea1);   
            node=calculate_node(x,y,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x,y,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x=insect1(1);y=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa3],[TX(2),bb3],'k');
               plot([aa3,aa1],[bb3,bb1],'k');
               plot([aa1,aa2],[bb1,bb2],'k');
               xx=aa2;yy=bb2;
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end    
 end 
 
 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
         for k=1:1

             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa4=x1;bb4=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 
            aa3=x1;bb3=y1;  

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 

            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa1=x1;bb1=y1;
            aa2=x;bb2=y;
            x2=x+cos(thea1);y2=y+sin(thea1);   
            node=calculate_node(x,y,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x,y,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x=insect1(1);y=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa4],[TX(2),bb4],'k');
               plot([aa4,aa3],[bb4,bb3],'k');
               plot([aa3,aa1],[bb3,bb1],'k');
               plot([aa1,aa2],[bb1,bb2],'k');
               xx=aa2;yy=bb2;
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end    
 end         


 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa5=x1;bb5=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa4=x1;bb4=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 
            aa3=x1;bb3=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 

            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa2=x1;bb2=y1;
            aa1=x;bb1=y;
            value1=calculate(x,y,RX(1),RX(2),building);
            if value1==1
                plot([TX(1),aa5],[TX(2),bb5],'r');
                plot([aa5,aa4],[bb5,bb4],'r');
                plot([aa4,aa3],[bb4,bb3],'r');
                plot([aa3,aa2],[bb3,bb2],'r');
                plot([aa2,aa1],[bb2,bb1],'r');
                plot([aa1,RX(1)],[bb1,RX(2)],'r');
            end
    end      
 end

（2）calculate.m詳細代碼

function value=calculate(x2,y2,x3,y3,build)
    N=17;
    k=1;mm=0;
    value=0;
    while k >=1
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            X1=[x3,y3];Y1=[x2,y2];X2=build{k}(p,:);Y2=build{k}(p+1,:);


            if X1(1)==Y1(1)
                X=X1(1);
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b2=X2(2)-k2*X2(1); 
                Y=k2*X+b2;
            end
            if X2(1)==Y2(1)
                X=X2(1);
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                b1=X1(2)-k1*X1(1);
                Y=k1*X+b1;
            end
            if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b1=X1(2)-k1*X1(1);
                b2=X2(2)-k2*X2(1);
                if k1==k2
                    continue
                else
                    X=(b2-b1)/(k1-k2);
                    Y=k1*X+b1;
                end
            end


            Q=([X,Y]-X1)*(Y1-[X,Y])';
            W=([X,Y]-X2)*(Y2-[X,Y])';
            if (Q >0.001)&&(W >0.001)
                mm=1;
            end
        end 
        if mm==1
            mm=0;
            k=1;
            break
        end    
        k=k+1;     
        if k >N
            break
        end
    end
    if k >N
        value=1;
        k=1;
    end  
   end

（3）calculate_node.m詳細代碼

function node=calculate_node(x1,y1,build)
    N=17;m=0;node=0;
    for k=1:N
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            x2=build{k}(p,1);y2=build{k}(p,2);
         if abs(x1-x2)<=0.1&&abs(y1-y2)<=0.1
             node=k;
             m=1;
             break
         end
         if m==1
             break 
         end
        end  
    end
end

（4）calculate_point.m詳細代碼

function [x1]=calculate_point(x2,y2,build)
    N=17;
    k=1;m=1;mm=0;dd=0;
    node=calculate_node(x2,y2,build);
    for i=1:N
        if i==node
            continue
        end
        row_column=size(build{i});
        row=row_column(1);

        value=findway(x2,y2,build);
        if value==i
            continue
        end

        for j=1:row-1
            x=build{i}(j,1);
            y=build{i}(j,2);
            while k >=1
                row_column1=size(build{k});
                row1=row_column1(1);
                for p=1:row1-1
                    X1=[x,y];Y1=[x2,y2];X2=build{k}(p,:);Y2=build{k}(p+1,:);

                    if (X1(1)==Y1(1)&&X1(2)==Y1(2))  ||  (X2(1)==Y2(1)&&X2(2)==Y2(2))
                        continue
                    end
                    if X1(1)==Y1(1)
                        X=X1(1);
                        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                        b2=X2(2)-k2*X2(1); 
                        Y=k2*X+b2;
                    end
                    if X2(1)==Y2(1)
                        X=X2(1);
                        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                        b1=X1(2)-k1*X1(1);
                        Y=k1*X+b1;
                    end
                    if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
                        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                        b1=X1(2)-k1*X1(1);
                        b2=X2(2)-k2*X2(1);
                        if k1==k2
                            continue
                        else
                            X=(b2-b1)/(k1-k2);
                            Y=k1*X+b1;
                        end
                    end

                    Q=([X,Y]-X1)*(Y1-[X,Y])';
                    W=([X,Y]-X2)*(Y2-[X,Y])';
                    if (Q >0.001)&&(W >0.001)
                        mm=1;
                    end
                end 
                if mm==1
                    mm=0;
                    k=1;
                    break
                end    
                k=k+1;     
                if k >N
                    break
                end
            end
            if k >N
                for ii=1:N
                    row_column11=size(build{ii});
                    row11=row_column11(1);
                    for pp=1:row11-1
                        x22=build{ii}(pp,1);y22=build{ii}(pp,2);
                        QQ=([x22,y22]-Y1)*([x,y]-[x22,y22])'/norm([x22,y22]-Y1)/norm([x,y]-[x22,y22]);
                        if QQ==1
                            dd=1;
                            break
                        end
                    end
                    if dd==1
                            break
                    end
                end
                if dd==0
                    x1(m,1)=x;x1(m,2)=y;
                    m=m+1;
                end
                k=1;
                dd=0;
            end          
        end
    end
end

（5）findway.m詳細代碼

function value=findway(xx,yy,build)
    N=17;
    value=0;
    for i=1:N
        row_column=size(build{i});
        row=row_column(1);
        for j=1:row-1
            x=build{i}(j,1);
            y=build{i}(j,2);
            if x==xx&&y==yy
                value=i;
                break
            end
        end
        if value==i
            break
        end
    end
end

（6）insection_point.m詳細代碼

function [point1,point2,point3,a,b,c,d]=insection_point(x2,y2,x3,y3,build)
    N=17;min_num=10^10;
    point1(1)=0;point1(2)=0;point2(1)=0;point2(2)=0;point3(1)=0;point3(2)=0;a=0;b=0;c=0;d=0;
    for k=1:N
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            X1=[x3,y3];Y1=[x2,y2];X2=build{k}(p,:);Y2=build{k}(p+1,:);
            if X1(1)==Y1(1)
                X=X1(1);
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b2=X2(2)-k2*X2(1); 
                Y=k2*X+b2;
            end
            if X2(1)==Y2(1)
                X=X2(1);
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                b1=X1(2)-k1*X1(1);
                Y=k1*X+b1;
            end
            if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b1=X1(2)-k1*X1(1);
                b2=X2(2)-k2*X2(1);
                if k1==k2
                    continue
                else
                    X=(b2-b1)/(k1-k2);
                    Y=k1*X+b1;
                end
            end


            Q=(X1-Y1)*([X,Y]-Y1)';
            W=([X,Y]-X2)*(Y2-[X,Y])';
            if (Q >0.001)&&(W >0.001)
                num1=(x2-X)^2+(y2-Y)^2;
               if min_num >=num1
                   min_num=num1;
                   point1=[X,Y];
                   a=build{k}(p,1);b=build{k}(p,2);
                   c=build{k}(p+1,1);d=build{k}(p+1,2);
                       if a==c
                            point2=[2*a-x2,y2];
                            point3=point2;
                            point2=2.*point1-point2;
                       end
                       if b==d
                            point2=[x2,2*d-y2];
                            point3=point2;
                            point2=2.*point1-point2;
                       end
               end              
            end            
        end  
    end
end

（7）line_to_point.m詳細代碼

function [x,y]=line_to_point(x1,y1,x2,y2,build)
    x=0;y=0;
    N=17;
    for k=1:N
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            x3=build{k}(p,1);y3=build{k}(p,2);
            value1=calculate(x1,y1,x3,y3,build);
            zhi=([x3,y3]-[x1,y1])*([x2,y2]-[x1,y1])'/norm([x3,y3]-[x1,y1])/norm([x2,y2]-[x1,y1]);
            if zhi >=0.99&&zhi<=1&&value1==1
                x=x3;y=y3;
            end
        end
    end
end

（8）point_in.m詳細代碼

function [value1,value2]=point_in(x2,y2,x3,y3,a,b,c,d)
    value1=0;value2=0;
    X1=[x3,y3];Y1=[x2,y2];X2=[a,b];Y2=[c,d];
    if X1(1)==Y1(1)
        X=X1(1);
        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
        b2=X2(2)-k2*X2(1); 
        Y=k2*X+b2;
    end
    if X2(1)==Y2(1)
        X=X2(1);
        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
        b1=X1(2)-k1*X1(1);
        Y=k1*X+b1;
    end
    if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
        b1=X1(2)-k1*X1(1);
        b2=X2(2)-k2*X2(1);
        if k1==k2
            value1=0;value2=0;
        else
            X=(b2-b1)/(k1-k2);
            Y=k1*X+b1;
        end
    end
    Q=([X,Y]-X1)*(Y1-[X,Y])';
    W=([X,Y]-X2)*(Y2-[X,Y])';
    if (Q >0.001)&&(W >0.001)
       value1=X;value2=Y;
    end
end