bspline_2.m

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

%************************************************************************ 
%* 文件名：    Spline_B.m
%* 函数名：    Spline_B
%* 参数：      Info,U,V（Info为原始轮廓信息，包括是否是轮廓以及点数量；U，V为X，Y坐标数值数组10*30）
%* 功能描述:   实现原始轮廓数据的二次B样条平滑，根据特征多边形拟合闭合曲线
%* 返回值：    Spline_Info,Spline_U,Spline_V（Spline_Info为拟合后的轮廓信息，包括是否是轮廓以及点数量；Spline_U，Spline_V为X，Y坐标数值数组10*3000）
%* 作 者：     王刚 · [e-mail] · 哈尔滨理工大学 · 研究生05-1班
%* 完成时间：  2006年05月30日
%* 版本号：    1.0
%************************************************************************ 
function [Spline_Info,Spline_U,Spline_V,Spline_R]=BSpline_2(Info,U,V,R);
%U，R为原来的X，Y点。。。
Spline_Info_R=zeros(10,2);
Spline_U_R=zeros(10,30);
Spline_V_R=zeros(10,30);
Spline_R_R=zeros(10,30);

Spline_Info_R=Info;
Spline_U_R=U;
Spline_V_R=V;                                                       
Spline_R_R=R;                                                       %Spline_Info_R,Spline_U_R,Spline_V_R记录Info,U,V的值，不用删除也可以


for loop_R=1:10                                                     %扫描整个平面的轮廓线从1到10
if (Info(loop_R,1)==1)                                              %若是轮廓线逻辑Info(loop_R,1)即为1，则处理，否则放弃
    Spline_U_R(loop_R,(Info(loop_R,2)+1))=U(loop_R,1);
    Spline_V_R(loop_R,(Info(loop_R,2)+1))=V(loop_R,1);              %数组末尾处插入第一点数据，也就是轮廓的第二点数据值
    Spline_R_R(loop_R,(Info(loop_R,2)+1))=R(loop_R,1);
    x_temp=Spline_U_R(loop_R,Info(loop_R,2));
    y_temp=Spline_V_R(loop_R,Info(loop_R,2));                       %存储临时变量记录原始数组最后一点的数值，也就是轮廓的第一点数据值
    z_temp=Spline_R_R(loop_R,Info(loop_R,2));
    Spline_Info_R(loop_R,2)=Spline_Info_R(loop_R,2)+2;              %为拟合B样条，增加两点数据，所以点信息要增加2
    
    for loop_Move=Spline_Info_R(loop_R,2):-1:2
        Spline_U_R(loop_R,loop_Move)=Spline_U_R(loop_R,loop_Move-1);
        Spline_V_R(loop_R,loop_Move)=Spline_V_R(loop_R,loop_Move-1);
        Spline_R_R(loop_R,loop_Move)=Spline_R_R(loop_R,loop_Move-1);
    end    
    Spline_U_R(loop_R,1)=x_temp;
    Spline_V_R(loop_R,1)=y_temp;                                     %数据从后向前移动，最后第一点数值（x_tempy_temp）插入
    Spline_R_R(loop_R,1)=z_temp;
end
end
%************************************************************************ 
%* 经过以上处理，已经计算出了B样条关键点矩阵Spline_U_R，Spline_V_R
%* 以下为算法实现核心部分
%* 算法描述:  P（t）=(t^2 t 1)*1/2*[1 -2 1;-2 2 0;1 1 0]*[p0 p1 p2],其中0<=t<=1
%* 算法实现:  根据Pi=(xi,yi)为二次B样条曲线特征多边形三顶点矢量，以上抽象公式表示成分量形式，即以X(t),Y(t)的形式出现
%* 作 者：    王刚 · [e-mail] · 哈尔滨理工大学 · 研究生05-1班 
%* 完成时间： 2006/05/28 PM 9：30
%************************************************************************
CONST_POINT=5                                                       %CONST_POINT为点数扩大100倍，若减少计算量，可以适当控制扩大倍数，比如20，50等
Spline_U=zeros(10,30*CONST_POINT);
Spline_V=zeros(10,30*CONST_POINT);
Spline_R=zeros(10,30*CONST_POINT);
Spline_Info=zeros(10,2);                                              %Spline_Info(k,1)记录是否是轮廓，Spline_Info(k,2)记录点的个数
for loop=1:10
if (Info(loop,1)==1)
    Spline_Info(loop,2)=Info(loop,2)*CONST_POINT;
    Spline_Info(loop,1)=1;                                            
end
end

spline_n=CONST_POINT;
spline_dt=1/spline_n;                                                  %微分量spline_dt

for loop=1:10
    loop_Sum=0;                                                        %定义计数器loop_Sum，计算完后增加点值应该为N*spline_n
    if Spline_Info_R(loop,1)==1
    
        for loop_R=1:Spline_Info_R(loop,2)-2
          spline_a0=(Spline_U_R(loop,loop_R)+Spline_U_R(loop,loop_R+1))/2;
          spline_a1=Spline_U_R(loop,loop_R+1)-Spline_U_R(loop,loop_R);
          spline_a2=(Spline_U_R(loop,loop_R)-2*Spline_U_R(loop,loop_R+1)+Spline_U_R(loop,loop_R+2))/2;

          spline_b0=(Spline_V_R(loop,loop_R)+Spline_V_R(loop,loop_R+1))/2;
          spline_b1=Spline_V_R(loop,loop_R+1)-Spline_V_R(loop,loop_R);
          spline_b2=(Spline_V_R(loop,loop_R)-2*Spline_V_R(loop,loop_R+1)+Spline_V_R(loop,loop_R+2))/2;
          
          spline_c0=(Spline_R_R(loop,loop_R)+Spline_R_R(loop,loop_R+1))/2;
          spline_c1=Spline_R_R(loop,loop_R+1)-Spline_R_R(loop,loop_R);
          spline_c2=(Spline_R_R(loop,loop_R)-2*Spline_R_R(loop,loop_R+1)+Spline_R_R(loop,loop_R+2))/2;
                                                                        %以上spline_a0...spline_b2分别为曲线参数计算                                                                        
          for loop_Inter=1:spline_n                                      
               spline_t=loop_Inter*spline_dt;
               loop_Sum=loop_Sum+1;
               Spline_U(loop,loop_Sum)=spline_a0+spline_a1*spline_t+spline_a2*spline_t*spline_t;
               Spline_V(loop,loop_Sum)=spline_b0+spline_b1*spline_t+spline_b2*spline_t*spline_t;
               Spline_R(loop,loop_Sum)=spline_c0+spline_c1*spline_t+spline_c2*spline_t*spline_t;
          end                                                           %以上计算总拟合数据点，(Spline_Info_R(loop,2)-2)*spline_n个
        end
%************************************************************************ 
%* 经过以上处理，已经计算出了B样条关键点矩阵Spline_U，Spline_V
%* 以下为算法图形显示部分,为实现方便，与以上点计算模块一同嵌套在LOOP循环内，但不符合结构化的设计风格
%* 实现原理:使用plot函数,根据数值点总数建立图象显示阵列New_Spline_U，New_Spline_V，需要把原来矩阵多余的0去掉，否则显示不正常
%* 作 者：    王刚 · [e-mail] · 哈尔滨理工大学 · 研究生05-1班 
%* 完成时间： 2006/05/29 AM 10：36
%************************************************************************ 
          New_Spline_U=zeros(1,Spline_Info(loop,2));
          New_Spline_V=zeros(1,Spline_Info(loop,2));                    %define the array of figure show
          New_Spline_R=zeros(1,Spline_Info(loop,2));
          
          for New_Array=1:Spline_Info(loop,2)
              New_Spline_U(1,New_Array)=Spline_U(loop,New_Array);
              New_Spline_V(1,New_Array)=Spline_V(loop,New_Array);
              New_Spline_R(1,New_Array)=Spline_R(loop,New_Array);
          end                                                           %compute
          
          plot3(New_Spline_U(1,:),New_Spline_R(1,:),New_Spline_V(1,:));                    %show
    end
end