%% 数据准备与参数初始化 clear; clc; dataPath = 'D:\xt\新光源-第一次试验\1\'; % 24张图像的存储路径 imageFiles = dir(fullfile(dataPath, '*.jpg')); % 根据实际格式修改 numLights = 24; imageSize = [3000, 4096, 3]; % 根据实际尺寸修改 % 光源参数配置 latitude_angles = [15, 55, 75]; % 低中高纬度倾角 longitude_steps = 0:45:315; % 经度划分 light_directions = zeros(numLights, 3); % 存储光源方向向量 %% 步骤1：计算光源方向向量 idx = 1; for lat = latitude_angles theta = deg2rad(lat); for lon = longitude_steps phi = deg2rad(lon); light_directions(idx,:) = [sin(theta)*cos(phi),... sin(theta)*sin(phi),... cos(theta)]; idx = idx + 1; end end %% 步骤2：加载图像数据并预处理 rawImages = zeros([imageSize, numLights], 'uint8'); for i = 1:numLights rawImages(:,:,:,i) = imread(fullfile(dataPath, imageFiles(i).name)); end % 转换为双精度并线性化（假设已做相机响应校准） linearImages = im2double(rawImages); %% 步骤3：分离漫反射与镜面反射（RPCA实现） % 注意：需要安装MATLAB的统计和机器学习工具箱 lambda = 1/sqrt(max(size(linearImages))); [r, c, ch, n] = size(linearImages); separatedResults = zeros(r, c, ch, 2, n); % 存储分离结果 for ch = 1:3 % 分别处理RGB通道 channelData = squeeze(linearImages(:,:,ch,:)); for x = 1:r for y = 1:c [L,S] = rpca(reshape(channelData(x,y,:), [], numLights), 'Lambda', lambda); separatedResults(x,y,ch,1,:) = L; separatedResults(x,y,ch,2,:) = S; end end end diffuseComponents = separatedResults(:,:,:,1,:); % 漫反射分量 specularComponents = separatedResults(:,:,:,2,:); % 镜面反射分量 %% 步骤4：光度立体法计算法向量与漫反射系数 L = light_directions'; % 光源方向矩阵（3x24） diffuseParams = zeros(r, c, 4); % 存储[p, q, k_d, confidence] for x = 1:r for y = 1:c % 提取当前像素的漫反射观测值（RGB均值） I_d = mean(squeeze(diffuseComponents(x,y,:,1,:)), 1)'; % 求解线性方程组：I_d = k_d * L' * n A = L'; b = I_d; x_hat = A \ b; % 最小二乘解 % 法向量归一化 n = x_hat / norm(x_hat); p = -n(1)/n(3); % 表面梯度p分量 q = -n(2)/n(3); % 表面梯度q分量 k_d = norm(x_hat); % 漫反射系数 % 存储结果 diffuseParams(x,y,:) = [p, q, k_d, 1/(cond(A'*A))]; end end %% 步骤5：镜面反射系数估计 m = 50; % 高光指数（需根据材质调整） specularCoeff = zeros(r, c); for x = 1:r for y = 1:c % 获取当前像素法向量 n = [-diffuseParams(x,y,1), -diffuseParams(x,y,2), 1]; n = n / norm(n); % 寻找最佳镜面反射光源 maxSpec = 0; for i = 1:numLights l = light_directions(i,:)'; v = [0;0;1]; % 视角方向 % 计算镜面反射方向 r = 2*dot(n,l)*n - l; specTerm = (dot(r,v))^m / max(dot(n,l), 1e-6); % 取RGB通道最大值 currentSpec = max(specularComponents(x,y,:,2,i)); if specTerm > 1e-6 candidate = currentSpec * dot(n,l) / specTerm; maxSpec = max(maxSpec, candidate); end end specularCoeff(x,y) = maxSpec; end end %% 步骤6：计算梯度场 gradientMagnitude = zeros(r, c); gradientDirection = zeros(r, c); for x = 1:r for y = 1:c p = diffuseParams(x,y,1); q = diffuseParams(x,y,2); gradientMagnitude(x,y) = sqrt(p^2 + q^2); gradientDirection(x,y) = atan2(q, p); end end %% 结果可视化 figure; subplot(231), imshow(mean(diffuseComponents(:,:,:,1,1),3)), title('漫反射分量'); subplot(232), imshow(specularCoeff, []), title('镜面反射系数'); subplot(233), imshow(gradientMagnitude, []), title('梯度幅值'); subplot(234), quiver(diffuseParams(:,:,2), diffuseParams(:,:,1)), title('梯度方向场'); subplot(235), imshow(diffuseParams(:,:,3), []), title('漫反射系数'); subplot(236), polarhistogram(gradientDirection(:), 20), title('梯度方向分布');这段程序在运行的时候会报错：函数或变量 'rpca' 无法识别。 出错 Untitled (line 45) [L,S] = rpca(reshape(channelData(x,y,:), [], numLights), 'Lambda', lambda); >>