% 批量点云数据处理与统计分析程序
clear all
close all
clc;
format short
% 定义要处理的文件列表
% 定义要处理的文件列表(按需求添加/删除)
files = {
'a3.asc',
'A7.asc',
'c3.asc',
'B16.asc',
'c3.asc',
'c9.asc',
'c10.asc',
'c32.asc',
'c92.asc',
'c102.asc',
'D2.asc',
'D4.asc',
'D5.asc',
'D6.asc',
'D10.asc',
'D13.asc'
'D15.asc',
'D42.asc',
'D52.asc',
'D102.asc',
'D132.asc'
}; % 可继续在下方按格式添加更多文件}; % 可以根据需要添加更多文件
% 检查文件是否存在
existing_files = {};
for i = 1:length(files)
if exist(files{i}, 'file')
existing_files{end+1} = files{i};
else
fprintf('警告: 文件 %s 不存在\n', files{i});
end
end
if isempty(existing_files)
error('没有找到任何有效的ASC文件');
end
fprintf('将处理以下文件:\n');
for i = 1:length(existing_files)
fprintf('%d. %s\n', i, existing_files{i});
end
% 初始化结果存储结构
file_names = {}; % 存储文件名
results_data = []; % 存储结果数据
for file_idx = 1:length(existing_files)
fprintf('\n开始处理文件: %s (%d/%d)\n', existing_files{file_idx}, file_idx, length(existing_files));
try
% 加载当前文件
[x,y,z]=textread(existing_files{file_idx},'%f%f%f');
catch ME
fprintf('加载文件 %s 失败: %s\n', existing_files{file_idx}, ME.message);
continue;
end
% 检查数据是否为空或不完整
if isempty(x) || isempty(y) || isempty(z) || length(x) ~= length(y) || length(x) ~= length(z)
fprintf('文件 %s 数据不完整或格式错误,跳过...\n', existing_files{file_idx});
continue;
end
% %画出点云图像
% figure('Name', ['点云三维图像 - ', existing_files{file_idx}]);
% pcshow([x,y,z])
% colorbar;
% title(['点云三维图像 - ', existing_files{file_idx}]);
% xlabel('X');
% ylabel('Y');
% zlabel('Z');
% %直接生成三维曲面图
% tri=delaunay(x,y);
% figure('Name', ['三维曲面图 - ', existing_files{file_idx}]);
% trimesh(tri,x,y,z);
% colorbar;
% xlabel('x');
% ylabel('y');
% zlabel('z');
% --------------------------------------
%最小二乘法求解点云基平面
U=[x y ones(length(x),1)];
R=(pinv(U'*U)*(U')*z);
%基平面的系数向量(z=Ax+By+C)
A=R(1);
B=R(2);
C=R(3);
%建立相对于基平面的坐标
xi=x;
yi=y;
zi=A*xi+B*yi+C;
% %画出基平面图
% figure('Name', ['基平面图 - ', existing_files{file_idx}]);
% pcshow([x,y,z])
% hold on
% plot3(xi,yi,zi)
% title(['基平面图 - ', existing_files{file_idx}]);
% xlabel('X');
% ylabel('Y');
% zlabel('Z');
% %旋转基平面为xy平面
M=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 -C 1];
AA=sqrt(B.^2+1);BB=sqrt(A.^2+B.^2+1);CC=sqrt(B.^2+1);
R=[AA./BB 0 A./BB 0;-(A.*B)./(AA.*BB) 1./CC B./BB 0;A./(AA.*BB) B./CC -1./BB 0;0 0 0 1];
P=[x y z ones(length(x),1)];
pu=P*M*R;
x=pu(:,1);y=pu(:,2);z=-pu(:,3);
%最小二乘法求解点云基平面
U=[x y ones(length(x),1)];
R=(pinv(U'*U)*(U')*z);
%基平面的系数向量(z=Ax+By+C)
A=R(1);
B=R(2);
C=R(3);
%建立相对于基平面的坐标
xi=x;
yi=y;
zi=A*xi+B*yi+C;
% %画出处理后点云图
% figure('Name', ['处理后点云图 - ', existing_files{file_idx}]);
% pcshow([x,y,z])
% colorbar;
% xlabel('X');
% ylabel('Y');
% zlabel('Z');
% %画出处理后基平面图
% figure('Name', ['处理后基平面图 - ', existing_files{file_idx}]);
% pcshow([x,y,z])
% colorbar;
% hold on
% plot3(xi,yi,zi)
% title(['处理后基平面图 - ', existing_files{file_idx}]);
% xlabel('X');
% ylabel('Y');
% zlabel('Z');
% zi=z-zi;
zi=z-zi;
%%-------------------------------------------
% xy等间距化,绘制等间距点云图
x0=min(x)+0.5;x00=max(x)-0.5;y0=min(y)+0.5;y00=max(y)-0.5;
grid_size=0.2;
[xi,yi]=meshgrid(x0:grid_size:x00,y0:grid_size:y00);
%起点和终点的选择很重要,保证zi不出现NaN
zi=griddata(x,y,z,xi,yi,'cubic');
% %画出插值后点云图
% figure('Name', ['插值后点云图 - ', existing_files{file_idx}]);
% pcshow([xi(:) yi(:) zi(:)])
% colorbar;
% xlabel('X');
% ylabel('Y');
% zlabel('Z');
% %画出表面图
% figure('Name', ['表面图 - ', existing_files{file_idx}]);
% surf(xi,yi,zi);
% colorbar;
% shading interp;
% 选择一个剖面进行二维分析
% 这里我们选择Y方向的一个剖面,例如Y=mean(yi(:))
y_mid = mean(yi(:));
[~, idx] = min(abs(yi(:) - y_mid));
% 选择这个Y值对应的X-Z剖面
X1 = xi(idx, :);
Z1 = zi(idx, :);
% 移除NaN值
valid_idx = ~(isnan(X1) | isnan(Z1));
X1 = X1(valid_idx);
Z1 = Z1(valid_idx);
% 检查是否有足够的有效数据点进行分析
if length(X1) < 2
fprintf('文件 %s 有效数据点不足,跳过二维分析...\n', existing_files{file_idx});
% 为该文件设置默认值
Z2 = NaN; Rp = NaN; SF = NaN; C_2D = NaN; ThetamaxC1_2D = NaN; max_angle_deg = NaN;
else
% 验证输入数据
if length(X1) ~= length(Z1) || length(X1) < 2
error('X1和Z1必须是相同长度的数组,且长度至少为2');
end
% 初始化统计参数
sum_dert = 0; % Z方向变化量平方与X方向变化量比值之和
sum_sqt = 0; % 线段长度之和
sum_SF = 0; % Z方向变化量平方与X方向变化量乘积之和
total_L = 0; % X方向总长度
theta_array = []; % 倾角数组(弧度)
% 计算基本统计参数
for i = 1:length(X1)-1
delta_z = Z1(i+1) - Z1(i);
delta_x = X1(i+1) - X1(i);
if delta_x == 0
warning('发现X方向增量为0,跳过该点');
continue;
end
dert = (delta_z)^2 / delta_x;
dertsqt = sqrt(delta_z^2 + delta_x^2);
dertSF = delta_z^2 * delta_x;
sum_dert = sum_dert + dert;
sum_sqt = sum_sqt + dertsqt;
sum_SF = sum_SF + dertSF;
total_L = total_L + delta_x;
% 计算倾角(仅当Z方向变化为正时)
if delta_z >= 0 && dertsqt > 0
theta = asin(delta_z / dertsqt);
theta_array = [theta_array, theta];
else
theta_array = [theta_array, 0];
end
end
% 计算基本参数
avg_slope = sum_dert / total_L;
Z2 = sqrt(avg_slope); % 平均梯度的平方根
Rp = sum_sqt / total_L; % 平均线段长度
SF = sum_SF / total_L; % 曲线形状因子
% 将倾角转换为度
theta_degrees = theta_array * 180 / pi;
% 计算曲线总长度
total_curve_length = 0;
for i = 1:length(X1)-1
segment_length = sqrt((Z1(i+1)-Z1(i))^2 + (X1(i+1)-X1(i))^2);
total_curve_length = total_curve_length + segment_length;
end
% 计算不同倾角阈值下的曲线长度
angle_thresholds = 0:2:25; % 从0到25度,间隔2度
angle_thresholds_rad = angle_thresholds * pi / 180; % 转换为弧度
length_above_threshold = zeros(size(angle_thresholds));
for j = 1:length(angle_thresholds_rad)
current_length = 0;
for i = 1:length(X1)-1
current_slope = (Z1(i+1) - Z1(i)) / (X1(i+1) - X1(i));
current_angle = atan(abs(current_slope)); % 使用绝对值角度
if current_angle > angle_thresholds_rad(j)
segment_length = sqrt((Z1(i+1)-Z1(i))^2 + (X1(i+1)-X1(i))^2);
current_length = current_length + segment_length;
end
end
length_above_threshold(j) = current_length;
end
% 计算最大倾角
slope_angles = zeros(1, length(X1)-1);
for i = 1:length(X1)-1
if (X1(i+1) - X1(i)) ~= 0
slope = (Z1(i+1) - Z1(i)) / (X1(i+1) - X1(i));
slope_angles(i) = atan(abs(slope)); % 使用绝对值角度
else
slope_angles(i) = pi/2; % 垂直线段
end
end
max_angle_rad = max(slope_angles);
max_angle_deg = max_angle_rad * 180 / pi;
% 归一化长度比
normalized_lengths = length_above_threshold / total_curve_length;
% 非线性拟合
L0 = normalized_lengths(1); % 倾角大于0度的曲线长度比例
% 定义拟合函数
modelfun = @(C,theta) (L0 * ((max_angle_deg - theta) / max_angle_deg).^C);
% 初始参数估计
beta0 = [1];
try
% 执行非线性拟合
fitted_params = nlinfit(angle_thresholds, normalized_lengths, modelfun, beta0);
C_2D = fitted_params;
% %绘制拟合图像
% figure('Name', ['二维拟合图 - ', existing_files{file_idx}]);
% plot(angle_thresholds, normalized_lengths, 'bo', fit_theta, fit_curve, 'r-', 'LineWidth', 2);
% xlabel('倾角 θ (度)');
% ylabel('L_θ/L_0');
% title(['曲线倾角分布拟合 - ', existing_files{file_idx}]);
% legend('原始数据点', '拟合曲线', 'Location', 'best');
% grid on;
catch ME
warning('文件 %s 拟合失败: %s', existing_files{file_idx}, ME.message);
C_2D = NaN;
end
% 计算最终参数
if ~isnan(C_2D)
ThetamaxC1_2D = max_angle_deg / (C_2D + 1);
else
ThetamaxC1_2D = NaN;
end
end
% 显示结果摘要
fprintf('%s二维曲线几何特征分析结果:\n', existing_files{file_idx});
fprintf('平均梯度平方根 (Z2): %.4f\n', Z2);
fprintf('平均线段长度 (Rp): %.4f\n', Rp);
fprintf('曲线形状因子 (SF): %.4f\n', SF);
fprintf('拟合参数 (C_2D): %.4f\n', C_2D);
fprintf('最大倾角: %.2f度\n', max_angle_deg);
fprintf('Thetamax/(C+1): %.4f\n', ThetamaxC1_2D);
% --------------------------------------------------------------------------------------------------------
%计算三维粗糙度统计参数视倾角thetas 求基平面的法向力n0(注意法向力的方向问题)
X=xi(:);Y=yi(:);Z=zi(:);
% 对插值网格进行三角剖分,但需要处理NaN值
valid_idx = ~isnan(Z);
X_valid = X(valid_idx);
Y_valid = Y(valid_idx);
Z_valid = Z(valid_idx);
% 检查是否有足够的有效点进行三角剖分
if length(X_valid) >= 3
tri_grid = delaunay(X_valid, Y_valid);
n0=[A,B,-1];
n0=n0/norm(n0);
sumI=0;
sumAA=0;
m=size(tri_grid,1);
for i=1:m
point1= [X_valid(tri_grid(i,1)) Y_valid(tri_grid(i,1)) Z_valid(tri_grid(i,1))];
point2= [X_valid(tri_grid(i,2)) Y_valid(tri_grid(i,2)) Z_valid(tri_grid(i,2))];
point3= [X_valid(tri_grid(i,3)) Y_valid(tri_grid(i,3)) Z_valid(tri_grid(i,3))];
P1=point1-point2;
P2=point1-point3;
n=cross(P1,P2);
% 修复N的计算:N = (n*n0')*n0,即n在n0方向上的投影
if norm(n) > 1e-10
N = (n * n0') * n0;
n1 = n - N; %n剪切面的投影向量
% 修复:检查A是否为0,避免除零错误
if abs(A) > 1e-10
S=[C./A 0 0];
s= S/norm(S);
else
% 当A接近0时,使用不同的方向向量
S=[1 0 0]; % 默认x方向
s= S/norm(S);
end
% 检查向量是否有效
if norm(s) > 1e-10 && norm(n1) > 1e-10
% 限制acos函数的输入范围,避免NaN
dot_product = dot(n0,n)/(norm(n0)*norm(n));
dot_product = max(-1, min(1, dot_product));
Beta=acos(dot_product);
dot_product2 = dot(s,n1)/(norm(s)*norm(n1));
dot_product2 = max(-1, min(1, dot_product2));
alfa=acos(dot_product2);
I=-tan(Beta).*cos(alfa);
% 检查I是否为有限值
if isfinite(I)
sumI=sumI+I.^2;
end
sumAA=sumAA+1/2*norm(n);
end
end
end
%单位度
if sumAA > 0
Thetas=atan(sqrt(sumI/m));
Thetas=Thetas*180/pi; % 转化为角度
else
Thetas = 0;
fprintf('警告:无法计算Thetas,sumAA为0\n');
end
else
fprintf('警告:有效点数不足,无法进行三角剖分\n');
Thetas = 0;
end
% -----------%%计算三维参数Z2s、Rs、Thetas---------------------------------
Nx=size(xi,1);
Ny=size(xi,2);
dertaxi=grid_size;
dertayi=grid_size;
sumZ2s=0;
sumAt=0;
%zi是一个1032*1021的矩阵,而Nx、Ny值为1032,所以一直提示超出矩阵维度,注意zi中含有NaN!!!;
for i= 1: Nx-1
for j= 1: Ny-1
% 检查是否为NaN
if ~isnan(zi(i+1,j+1)) && ~isnan(zi(i,j+1)) && ~isnan(zi(i+1,j)) && ~isnan(zi(i,j)) && ...
~isnan(zi(i+1,j+1)) && ~isnan(zi(i+1,j)) && ~isnan(zi(i,j+1)) && ~isnan(zi(i,j))
zi_row1=(zi(i+1,j+1 )-zi(i,j+1 )).^2;
zi_row2=(zi(i+1,j)-zi(i,j)).^2;
zi_col1=(zi(i+1,j+1)-zi(i+1,j)).^2;
zi_col2=(zi(i,j+1)-zi(i,j)).^2;
Sz2s=(zi_row1+zi_row2)./(2*dertaxi.^2)+(zi_col1+zi_col2)./(2*dertayi.^2);
Ats=sqrt(1+zi_row2./(dertaxi.^2)+zi_col2./(dertayi.^2));
sumZ2s=sumZ2s+Sz2s;
sumAt=sumAt+Ats;
end
end
end
if (Nx-1)*(Ny-1) > 0
Z2s=(sumZ2s./((Nx-1).*(Ny-1))).^(1/2);
else
Z2s = 0;
fprintf('警告:无法计算Z2s,分母为0\n');
end
At=dertaxi*dertayi*sumAt;
An=(Nx-1)*(Ny-1)*dertaxi*dertayi;
if An > 0
Rs=At/An;
else
Rs = 0;
fprintf('警告:无法计算Rs,An为0\n');
end
%-----------------------注释7.19-------------
%计算三维情况下的最后一个参数
% 先求曲面初始总面积
sumA0=At;
% 检查sumA0是否有效
if sumA0 == 0 || ~isfinite(sumA0)
fprintf('警告:sumA0无效,值为%.6f\n', sumA0);
%return; % 如果sumA0无效,直接返回
else
% 求大于不同倾角的曲线长度
belta=0:2*pi/180:40*pi/180; % 设置不同间距下的倾角---% theta的取值应小于最大倾角,否则拟合的C为虚数
% 重新使用插值网格的三角剖分
if length(X_valid) >= 3
m=size(tri_grid,1);
for j=1:length(belta)
sumA(j)=0;
for i=1:size(tri_grid,1)
point1= [X_valid(tri_grid(i,1)) Y_valid(tri_grid(i,1)) Z_valid(tri_grid(i,1))];
point2= [X_valid(tri_grid(i,2)) Y_valid(tri_grid(i,2)) Z_valid(tri_grid(i,2))];
point3= [X_valid(tri_grid(i,3)) Y_valid(tri_grid(i,3)) Z_valid(tri_grid(i,3))];
P1=point1-point2;
P2=point1-point3;
n=cross(P1,P2);
if norm(n) > 1e-10 % 检查法向量是否为零向量
% 计算N = (n*n0')*n0,即n在n0方向上的投影
N = (n * n0') * n0;
n1 = n - N; %n剪切面的投影向量
% 修复:检查A是否为0,避免除零错误
if abs(A) > 1e-10
S=[C/A 0 0];
s= S/norm(S);
else
% 当A接近0时,使用不同的方向向量
S=[1 0 0]; % 默认x方向
s= S/norm(S);
end
if norm(s) > 1e-10 && norm(n1) > 1e-10 % 检查s和n1是否为有效向量
% 限制acos函数的输入范围,避免NaN
dot_product = dot(n0,n)/(norm(n0)*norm(n));
dot_product = max(-1, min(1, dot_product));
Beta=acos(dot_product);
dot_product2 = dot(s,n1)/(norm(s)*norm(n1));
dot_product2 = max(-1, min(1, dot_product2));
alfa=acos(dot_product2);
I=-tan(Beta).*cos(alfa);
% 检查I是否为有限值
if isfinite(I) && I > belta(j) %theta大于 belta(j)
As0=1/2*norm(n);
sumA(j)=sumA(j)+As0;
else
sumA(j)=sumA(j);
end
else
sumA(j)=sumA(j);
end
else
sumA(j)=sumA(j);
end
end
end
sumA(end);
% 求最大thetam
thetam3D = zeros(size(tri_grid,1), 1);
for i=1:size(tri_grid,1)
point1= [X_valid(tri_grid(i,1)) Y_valid(tri_grid(i,1)) Z_valid(tri_grid(i,1))];
point2= [X_valid(tri_grid(i,2)) Y_valid(tri_grid(i,2)) Z_valid(tri_grid(i,2))];
point3= [X_valid(tri_grid(i,3)) Y_valid(tri_grid(i,3)) Z_valid(tri_grid(i,3))];
P1=point1-point2;
P2=point1-point3;
n=cross(P1,P2);
if norm(n) > 1e-10 % 检查法向量是否为零向量
% 计算N = (n*n0')*n0,即n在n0方向上的投影
N = (n * n0') * n0;
n1 = n - N; %n剪切面的投影向量
% 修复:检查A是否为0,避免除零错误
if abs(A) > 1e-10
S=[C/A 0 0];
s= S/norm(S);
else
% 当A接近0时,使用不同的方向向量
S=[1 0 0]; % 默认x方向
s= S/norm(S);
end
if norm(s) > 1e-10 && norm(n1) > 1e-10 % 检查s和n1是否为有效向量
% 限制acos函数的输入范围,避免NaN
dot_product = dot(n0,n)/(norm(n0)*norm(n));
dot_product = max(-1, min(1, dot_product));
Beta=acos(dot_product);
dot_product2 = dot(s,n1)/(norm(s)*norm(n1));
dot_product2 = max(-1, min(1, dot_product2));
alfa=acos(dot_product2);
I=-tan(Beta).*cos(alfa);
% 检查I是否为有限值
if isfinite(I)
thetam3D(i)=atan(abs(I)); % 使用绝对值避免复数
else
thetam3D(i) = 0; % 默认值
end
else
thetam3D(i) = 0; % 默认值
end
else
thetam3D(i) = 0; % 默认值
end
end
thetamax3D=max(thetam3D);
thetamax3D=thetamax3D*180/pi;% 转化为角度
% 开始拟合出C值
Atheta=sumA/sumA0;
theta=belta;
theta=theta*180/pi;% 转化为角度
A0=sumA(1)/sumA0;% sumL(1)是倾角大于0度的面积
modelfun = @(C,theta)(A0*((thetamax3D-theta)./thetamax3D).^C);% theta的取值应小于最大倾角,否则拟合的C为虚数
beta0 = [1];
beta = nlinfit(theta,Atheta,modelfun,beta0);
C_3D= beta;
% %绘制拟合图像
% figure('Name', ['三维拟合图 - ', existing_files{file_idx}]);
% plot(theta,Atheta,'^','Markersize',10,'MarkerFaceColor','b')
% hold on
% plot(x22,y22)
% xlabel('θ');
% ylabel('Aθ');
% legend('原始数据点' , '拟合后的曲线')
%计算thetamax./(C+1)的值
ThetamaxC1_3D=thetamax3D/(C_3D+1);
% %单位:度 ThetamaxC1_3D=ThetamaxC1_3D*180/pi
else
fprintf('无法进行插值网格的三角剖分,有效点数不足\n');
C_3D = NaN;
ThetamaxC1_3D = NaN;
end
end
% 保存当前文件的计算结果到数组中
current_results = [Z2, Rp, SF, C_2D, ThetamaxC1_2D, max_angle_deg, C_3D, ThetamaxC1_3D, Rs, Z2s];
file_names{end+1} = existing_files{file_idx}; % 存储文件名
results_data = [results_data; current_results]; % 存储结果数据
% 显示当前文件的最终结果
fprintf('\n%s最终统计分析结果:\n', existing_files{file_idx});
fprintf('Z2: %.6f\n', Z2);
fprintf('Rp: %.6f\n', Rp);
fprintf('SF: %.6f\n', SF);
fprintf('C_2D: %.6f\n', C_2D);
fprintf('ThetamaxC1_2D: %.6f\n', ThetamaxC1_2D);
fprintf('max_angle_deg: %.6f\n', max_angle_deg);
fprintf('C_3D: %.6f\n', C_3D);
fprintf('ThetamaxC1_3D: %.6f\n', ThetamaxC1_3D);
fprintf('Rs: %.6f\n', Rs);
fprintf('Z2s: %.6f\n', Z2s);
fprintf('----------------------------------------\n');
end
% 创建汇总结果的Excel表格
fprintf('\n正在生成汇总结果Excel表格...\n');
if ~isempty(results_data) && ~isempty(file_names)
% 确保file_names是列向量
if isrow(file_names)
file_names = file_names';
end
% 确保results_data的行数与file_names一致
if length(file_names) == size(results_data, 1)
% 使用table函数一次性创建表格
T = table(file_names, results_data(:,1), results_data(:,2), results_data(:,3), ...
results_data(:,4), results_data(:,5), results_data(:,6), ...
results_data(:,7), results_data(:,8), results_data(:,9), results_data(:,10), ...
'VariableNames', {'File', 'Z2', 'Rp', 'SF', 'C_2D', 'ThetamaxC1_2D', ...
'max_angle_deg', 'C_3D', 'ThetamaxC1_3D', 'Rs', 'Z2s'});
fprintf('汇总参数表格:\n');
disp(T);
% 保存到Excel文件
try
filename = '汇总统计分析结果.xlsx';
writetable(T, filename);
fprintf('汇总结果已保存到文件: %s\n', filename);
catch ME
fprintf('保存Excel文件失败: %s\n', ME.message);
fprintf('尝试保存为CSV格式...\n');
try
filename_csv = '汇总统计分析结果.csv';
writetable(T, filename_csv);
fprintf('汇总结果已保存到文件: %s\n', filename_csv);
catch ME2
fprintf('保存CSV文件也失败: %s\n', ME2.message);
end
end
else
fprintf('错误: file_names长度(%d)与results_data行数(%d)不匹配\n', ...
length(file_names), size(results_data, 1));
end
else
fprintf('没有有效的结果数据可以保存\n');
end
fprintf('\n所有文件处理完成!共处理 %d 个文件\n', length(existing_files));
% 显示处理完成信息
if ~isempty(results_data)
fprintf('\n处理完成的文件列表:\n');
for i = 1:length(file_names)
fprintf('%d. %s\n', i, file_names{i});
end
end
% 如果需要,可以显示结果表格
if ~isempty(results_data)
fprintf('\n最终结果表格预览:\n');
if size(results_data, 1) <= 10
disp(T);
else
disp(T(1:10, :)); % 只显示前10行
fprintf('... (还有 %d 行结果)\n', size(results_data, 1)-10);
end
end
这是我的岩石结构面粗糙度分析代码,检查有没有错误,并解释每段代码