gen_mose_bd_2.m

该代码提取图片某一行或者某一列信息

%180--360_MOSE
index_x15_MOSE_negative=find(PD_x15_MOSE(:,3)<0);
PD_x15_MOSE_negative=PD_x15_MOSE(index_x15_MOSE_negative,1:4);
[PD_x15_MOSE_negative_z_sort,K_MOSE]=sort(PD_x15_MOSE_negative(:,2),'descend');
theta_x15_MOSE_negative_sort=zeros(size(K_MOSE));
D_x15_MOSE_negative_z_sort=zeros(size(K_MOSE));
for ii=1:size(K_MOSE)
    D_x15_MOSE_negative_z_sort(ii)=PD_x15_MOSE_negative(K_MOSE(ii),4);
    theta_x15_MOSE_negative_sort(ii)=135+acos(PD_x15_MOSE_negative(K_MOSE(ii),2)/10)*180/pi;
end

theta_x15_MOSE_sort=[theta_x15_MOSE_positive_sort;theta_x15_MOSE_negative_sort];
D_x15_MOSE_sort=[D_x15_MOSE_positive_sort; D_x15_MOSE_negative_z_sort];
thetayD_MOSE=[theta_x15_MOSE_sort,D_x15_MOSE_sort];
[thetay_MOSE_sort, order_MOSE]=sort(thetayD_MOSE(:,1), 'ascend');
thetayD_MOSE_new=zeros(size(thetayD_MOSE(:,1)),2);
for hh=1:size(thetayD_MOSE(:,1))
    thetayD_MOSE_new(hh,:)=thetayD_MOSE(order_MOSE(hh),:);
end

max_MOSE_x15=max(D_x15_MOSE_sort)

ratio=max_FEM_x15/max_MOSE_x15


size(PD_x15_MOSE(:,3));
MOSE_xx=1:size(PD_x15_MOSE(:,3));
%Axis([0,360,0,1]); 
hold on, plot(thetayD_MOSE_new(:,1),thetayD_MOSE_new(:,2), 'r:','LineWidth',2);
MOSE_curve=zeros(size(PD_x15_MOSE(:,3)),2);
MOSE_curve=[thetayD_MOSE_new(:,1),thetayD_MOSE_new(:,2)]

%calculate the relative error between AMGM and MOSE on circle x=13
temp_DcloseMOSE=100;
temp_AbsError=0;
temp_FEMEnergy=0;
for i=1:size(FEM_curve(:,2))
    for j=1:size(MOSE_curve(:,2))
        if(abs(FEM_curve(i,1)-MOSE_curve(j,1)<=2))
            if(abs(FEM_curve(i,2)-MOSE_curve(j,2))<abs(FEM_curve(i,2)-temp_DcloseMOSE))
                temp_DcloseMOSE=MOSE_curve(j,2);
            end
        end        
    end
    temp_AbsError = temp_AbsError + (temp_DcloseMOSE-FEM_curve(i,2))^2;
    temp_FEMEnergy = temp_FEMEnergy + FEM_curve(i,2)^2;
end

RelativeError=temp_AbsError/temp_FEMEnergy
%relative error over

%calculate the relative error between AMGM and MOSE on the cylinder surface
vv=zeros(size(Dense_BP_Surface(:,1)),1);
temp_AbsErrorTotal=0;
temp_FEMEnergyTotal=0;
for i=1:size(Dense_BP_Surface(:,2))  
    tp_M = MOSE_data(:,1:3) - ones(size(MOSE_data(:,2)),1)*Dense_BP_Surface(i,1:3);   
    [vv(i),n] = min(sum((tp_M.').^2));   
    temp_AbsErrorTotal = temp_AbsErrorTotal + (Dense_BP_Surface(i,4)-MOSE_data(n,4))^2;
    temp_FEMEnergyTotal= temp_FEMEnergyTotal+ Dense_BP_Surface(i,4)^2;
end
RelativeErrorTotal=temp_AbsErrorTotal/temp_FEMEnergyTotal
% figure (2);
% plot(vv);grid on
% title('误差曲线');
%relative error over

%other ragions for rendering
% fid_left_lung = fopen('D:\KASKADE\KeVin\PS_MG\Multigrid_PS\AMGM\1S_OFF\llung.off','r');
% file_ext = fscanf(fid_left_lung,'%s',[1]);
% num_data_left_lung = fscanf(fid_left_lung,'%d',[3,1]);
% num_P_left_lung=num_data_left_lung(1)-1;
% triangle_num_left_lung = num_data_left_lung(2);
% dummy = fscanf(fid_left_lung,'%f',[3,1]);
% off_point_left_lung = fscanf(fid_left_lung,'%f',[3,num_P_left_lung]);
% off_point_left_lung=off_point_left_lung.';
% off_triangle0_left_lung = fscanf(fid_left_lung,'%f',[4,triangle_num_left_lung]);
% off_triangle0_left_lung=off_triangle0_left_lung.';
% off_triangle_left_lung=off_triangle0_left_lung(:,2:4);
% fclose(fid_left_lung);
% 
% 
% fid_right_lung = fopen('D:\KASKADE\KeVin\PS_MG\Multigrid_PS\AMGM\1S_OFF\rlung.off','r');
% file_ext = fscanf(fid_right_lung,'%s',[1]);
% num_data_right_lung = fscanf(fid_right_lung,'%d',[3,1]);
% num_P_right_lung=num_data_right_lung(1)-1;
% triangle_num_right_lung = num_data_right_lung(2);
% dummy = fscanf(fid_right_lung,'%f',[3,1]);
% off_point_right_lung = fscanf(fid_right_lung,'%f',[3,num_P_right_lung]);
% off_point_right_lung=off_point_right_lung.';
% off_triangle0_right_lung = fscanf(fid_right_lung,'%f',[4,triangle_num_right_lung]);
% off_triangle0_right_lung=off_triangle0_right_lung.';
% off_triangle_right_lung=off_triangle0_right_lung(:,2:4);
% fclose(fid_right_lung);
% 
% fid_heart = fopen('D:\KASKADE\KeVin\PS_MG\Multigrid_PS\AMGM\1S_OFF\heart.off','r');
% file_ext = fscanf(fid_heart,'%s',[1]);
% num_data_heart = fscanf(fid_heart,'%d',[3,1]);
% num_P_heart=num_data_heart(1)-1;
% triangle_num_heart = num_data_heart(2);
% dummy = fscanf(fid_heart,'%f',[3,1]);
% off_point_heart = fscanf(fid_heart,'%f',[3,num_P_heart]);
% off_point_heart=off_point_heart.';
% off_triangle0_heart = fscanf(fid_heart,'%f',[4,triangle_num_heart]);
% off_triangle0_heart=off_triangle0_heart.';
% off_triangle_heart=off_triangle0_heart(:,2:4);
% fclose(fid_heart);
% 
% fid_liver = fopen('D:\KASKADE\KeVin\PS_MG\Multigrid_PS\AMGM\1S_OFF\liver.off','r');
% file_ext = fscanf(fid_liver,'%s',[1]);
% num_data_liver = fscanf(fid_liver,'%d',[3,1]);
% num_P_liver=num_data_liver(1)-1;
% triangle_num_liver = num_data_liver(2);
% dummy = fscanf(fid_liver,'%f',[3,1]);
% off_point_liver = fscanf(fid_liver,'%f',[3,num_P_liver]);
% off_point_liver=off_point_liver.';
% off_triangle0_liver = fscanf(fid_liver,'%f',[4,triangle_num_liver]);
% off_triangle0_liver=off_triangle0_liver.';
% off_triangle_liver=off_triangle0_liver(:,2:4);
% fclose(fid_liver);
% 
% fid_bone = fopen('D:\KASKADE\KeVin\PS_MG\Multigrid_PS\AMGM\1S_OFF\bone.off','r');
% file_ext = fscanf(fid_bone,'%s',[1]);
% num_data_bone = fscanf(fid_bone,'%d',[3,1]);
% num_P_bone=num_data_bone(1)-1;
% triangle_num_bone = num_data_bone(2);
% dummy = fscanf(fid_liver,'%f',[3,1]);
% off_point_bone = fscanf(fid_bone,'%f',[3,num_P_bone]);
% off_point_bone=off_point_bone.';
% off_triangle0_bone = fscanf(fid_bone,'%f',[4,triangle_num_bone]);
% off_triangle0_bone=off_triangle0_bone.';
% off_triangle_bone=off_triangle0_bone(:,2:4);
% fclose(fid_bone);

% fid_source = fopen('D:\KASKADE\KeVin\PS_MG\Multigrid_PS\AMGM\1S_OFF\source.off','r');
% file_ext = fscanf(fid_source,'%s',[1]);
% num_data_source = fscanf(fid_source,'%d',[3,1]);
% num_P_source=num_data_source(1)-1;
% triangle_num_source = num_data_source(2);
% dummy = fscanf(fid_liver,'%f',[3,1]);
% off_point_source = fscanf(fid_source,'%f',[3,num_P_source]);
% off_point_source=off_point_source.';
% off_triangle0_source = fscanf(fid_source,'%f',[4,triangle_num_source]);
% off_triangle0_source=off_triangle0_source.';
% off_triangle_source=off_triangle0_source(:,2:4);
% fclose(fid_source);
%%over

% Dense_BP=zeros(num_P,1);
% for ii=1:num_P   
%     temp_sub = points_MOSE - ones(length(dense_MOSE),1)*off_point(ii,1:3);
%     [vv(ii),n] = min(sum((temp_sub.').^2));
%     if( min(sum((temp_sub.').^2))<0.001 )
%         Dense_BP(ii) =dense_MOSE(n)*ratio;
%     end
% end
% 
% %%在各部分的坐标值后面加上能量密度
% point_tecplot=zeros(num_P,4);
% point_tecplot(:,1:3)=off_point;
% point_tecplot(:,4)=Dense_BP;
% 
% point_tecplot_left_lung=zeros(num_P_left_lung,4);
% point_tecplot_left_lung(:,1:3)=off_point_left_lung;
% 
% point_tecplot_right_lung=zeros(num_P_right_lung,4);
% point_tecplot_right_lung(:,1:3)=off_point_right_lung;
% 
% point_tecplot_heart=zeros(num_P_heart,4);
% point_tecplot_heart(:,1:3)=off_point_heart;
% 
% point_tecplot_liver=zeros(num_P_liver,4);
% point_tecplot_liver(:,1:3)=off_point_liver;
% 
% point_tecplot_bone=zeros(num_P_bone,4);
% point_tecplot_bone(:,1:3)=off_point_bone;
% 
% % point_tecplot_source=zeros(num_P_source,4);
% % point_tecplot_source(:,1:3)=off_point_source;
% %%over
% 
% fid = fopen('D:\KASKADE\KeVin\PS_MG\Multigrid_PS\AMGM\2S\Boundary_MOSE_2S.dat','w');
%     % ============ 文件头 ================== %
% 	fprintf(fid,'%s\n','Boundary_MOSE ');
% 	fprintf(fid,'%s\n','VARIABLES = "X", "Y", "Z", "Density"');
%     
%     fprintf(fid,'%s','ZONE T="TRIANGLES", N=');
%     fprintf(fid,'%d',num_P);
%     fprintf(fid,'%s',', E=');
%     fprintf(fid,'%d',triangle_num);
%     fprintf(fid,'%s\n',', DATAPACKING=POINT, ZONETYPE=FETRIANGLE');
% 	fprintf(fid,'%9.6f\t %9.6f\t %9.6f\t %30.24f\n',point_tecplot.');             % 点坐标和能量值
%  	fprintf(fid,'\n \n');
% 	fprintf(fid,' %d\t %d\t %d\n',off_triangle.');          % 单元信息 
%     fprintf(fid,'\n \n');
%     
%     fprintf(fid,'%s','ZONE T="TRIANGLES", N=');
%     fprintf(fid,'%d',num_P_left_lung);
%     fprintf(fid,'%s',', E=');
%     fprintf(fid,'%d',triangle_num_left_lung);
%     fprintf(fid,'%s\n',', DATAPACKING=POINT, ZONETYPE=FETRIANGLE');
% 	fprintf(fid,'%9.6f\t %9.6f\t %9.6f\t %30.24f\n',point_tecplot_left_lung.');             % 点坐标和能量值
%  	fprintf(fid,'\n \n');
% 	fprintf(fid,' %d\t %d\t %d\n',off_triangle_left_lung.');          % 单元信息 
%     fprintf(fid,'\n \n');
%     
%     fprintf(fid,'%s','ZONE T="TRIANGLES", N=');
%     fprintf(fid,'%d',num_P_right_lung);
%     fprintf(fid,'%s',', E=');
%     fprintf(fid,'%d',triangle_num_right_lung);
%     fprintf(fid,'%s\n',', DATAPACKING=POINT, ZONETYPE=FETRIANGLE');
% 	fprintf(fid,'%9.6f\t %9.6f\t %9.6f\t %30.24f\n',point_tecplot_right_lung.');             % 点坐标和能量值
%  	fprintf(fid,'\n \n');
% 	fprintf(fid,' %d\t %d\t %d\n',off_triangle_right_lung.');          % 单元信息 
%     fprintf(fid,'\n \n');
%     
%     fprintf(fid,'%s','ZONE T="TRIANGLES", N=');
%     fprintf(fid,'%d',num_P_heart);
%     fprintf(fid,'%s',', E=');
%     fprintf(fid,'%d',triangle_num_heart);
%     fprintf(fid,'%s\n',', DATAPACKING=POINT, ZONETYPE=FETRIANGLE');
% 	fprintf(fid,'%9.6f\t %9.6f\t %9.6f\t %30.24f\n',point_tecplot_heart.');             % 点坐标和能量值
%  	fprintf(fid,'\n \n');
% 	fprintf(fid,' %d\t %d\t %d\n',off_triangle_heart.');          % 单元信息 
%     fprintf(fid,'\n \n');
%     
%     fprintf(fid,'%s','ZONE T="TRIANGLES", N=');
%     fprintf(fid,'%d',num_P_liver);
%     fprintf(fid,'%s',', E=');
%     fprintf(fid,'%d',triangle_num_liver);
%     fprintf(fid,'%s\n',', DATAPACKING=POINT, ZONETYPE=FETRIANGLE');
% 	fprintf(fid,'%9.6f\t %9.6f\t %9.6f\t %30.24f\n',point_tecplot_liver.');             % 点坐标和能量值
%  	fprintf(fid,'\n \n');
% 	fprintf(fid,' %d\t %d\t %d\n',off_triangle_liver.');          % 单元信息 
%     fprintf(fid,'\n \n');
%     
%     fprintf(fid,'%s','ZONE T="TRIANGLES", N=');
%     fprintf(fid,'%d',num_P_bone);
%     fprintf(fid,'%s',', E=');
%     fprintf(fid,'%d',triangle_num_bone);
%     fprintf(fid,'%s\n',', DATAPACKING=POINT, ZONETYPE=FETRIANGLE');
% 	fprintf(fid,'%9.6f\t %9.6f\t %9.6f\t %30.24f\n',point_tecplot_bone.');             % 点坐标和能量值
%  	fprintf(fid,'\n \n');
% 	fprintf(fid,' %d\t %d\t %d\n',off_triangle_bone.');          % 单元信息 
%     fprintf(fid,'\n \n');

%     fprintf(fid,'%s','ZONE T="TRIANGLES", N=');
%     fprintf(fid,'%d',num_P_source);
%     fprintf(fid,'%s',', E=');
%     fprintf(fid,'%d',triangle_num_source);
%     fprintf(fid,'%s\n',', DATAPACKING=POINT, ZONETYPE=FETRIANGLE');
% 	fprintf(fid,'%9.6f\t %9.6f\t %9.6f\t %30.24f\n',point_tecplot_source.');             % 点坐标和能量值
%  	fprintf(fid,'\n \n');
% 	fprintf(fid,' %d\t %d\t %d\n',off_triangle_source.');          % 单元信息 
%     fprintf(fid,'\n \n');
    
%fclose(fid);


% fid = fopen('D:\KASKADE\KeVin\PointS\Gen_Boundary\y8_22549\Boundary_sub.dat','w');
%     % ============ 文件头 ================== %
% 	fprintf(fid,'%s\n','Boundary_sub ');
% 	fprintf(fid,'%s\n','VARIABLES = "X", "Y", "Z", "Density"');
%     
%     fprintf(fid,'%s','ZONE T="TRIANGLES", N=');
%     fprintf(fid,'%d',num_P);
%     fprintf(fid,'%s',', E=');
%     fprintf(fid,'%d',triangle_num);
%     fprintf(fid,'%s\n',', DATAPACKING=POINT, ZONETYPE=FETRIANGLE');
% 	fprintf(fid,'%9.6f\t %9.6f\t %9.6f\t %30.24f\n',point_tecplot.');             % 点坐标和能量值
%  	fprintf(fid,'\n \n');
% 	fprintf(fid,' %d\t %d\t %d\n',off_triangle.');          % 单元信息 
%     fprintf(fid,'\n \n');
%     
% fclose(fid);%RelativeError1=,RelativeError2,RelativeError3=0.498,Relative
% Error4=.