【光学】基于非凸变分最小化优化经典光度立体附Matlab代码

 ​✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知，求助可私信。

🔥 内容介绍

光度立体 (Photometric Stereo) 作为一种无需依赖深度传感器即可恢复物体三维形状的计算机视觉技术，凭借其简单、高效的特性，在逆向工程、机器人视觉及文物数字化等领域得到了广泛应用。经典光度立体方法通常基于朗伯反射模型，通过分析同一场景下不同光照条件下的图像信息来重建物体的表面法向量，进而积分得到三维形状。然而，经典方法通常面临着非凸优化问题，导致解的唯一性难以保证，且对噪声较为敏感。本文将深入探讨基于非凸变分最小化优化策略来解决经典光度立体重建中存在的挑战。

经典光度立体方法的核心在于建立图像亮度与表面法向量之间的关系。假设物体表面满足朗伯反射模型，则像素点的亮度值与其表面法向量和光源方向的内积成正比。对于同一场景下M张不同光照条件下的图像，我们可以建立如下的线性方程组：

I = N L

其中，I 是一个 M × 1 的向量，表示像素点在M张图像中的亮度值；N 是一个 3 × 1 的向量，表示像素点的表面法向量；L 是一个 M × 3 的矩阵，其每一行表示一个光源的方向向量。

理想情况下，通过求解上述线性方程组即可得到表面法向量N。然而，由于实际采集图像中存在噪声以及朗伯反射模型的简化假设，直接求解该方程组往往会得到不准确甚至错误的结果。此外，该方程组是一个欠定方程组 (M < 3 时)，或者是一个超定方程组 (M > 3 时)，直接求解会导致解的不唯一性或者对噪声非常敏感。

为了解决这些问题，研究者们提出了各种优化策略。传统的最小二乘法是最常用的方法之一，其目标函数为：

min ||I - NL||^2

然而，最小二乘法容易受到噪声的影响，且其解的唯一性也难以保证。更重要的是，该目标函数是一个非凸函数，局部极小值的存在会严重影响优化结果的精度。

基于此，非凸变分最小化优化策略应运而生。其核心思想是通过引入正则化项来约束解的空间，从而避免陷入局部极小值，并提高对噪声的鲁棒性。常用的正则化项包括：

• L1正则化: L1正则化项可以有效地稀疏化解，去除冗余信息，提高对噪声的鲁棒性。其目标函数可以表示为

min ||I - NL||^2 + λ ||N||_1

其中，λ 是一个正则化参数，控制正则化项的权重。

• TV正则化: 全变分 (Total Variation, TV) 正则化项可以保持解的平滑性，抑制噪声的影响，尤其适用于表面光滑的物体。其目标函数可以表示为：

min ||I - NL||^2 + λ ||∇N||_1

其中，∇N 表示表面法向量的梯度。

• 高阶正则化: 为了更好地保持解的细节信息，可以使用高阶正则化项，例如二阶TV正则化。

这些非凸变分最小化问题通常采用迭代优化算法求解，例如梯度下降法、近端梯度法等。这些算法的收敛性及全局最优性的证明往往是极其困难的，需要结合具体的算法和正则化项进行分析。

此外，针对朗伯反射模型的局限性，研究者们也提出了许多改进方法，例如考虑镜面反射、漫反射与镜面反射混合等情况。这些改进方法通常会引入更加复杂的模型，使得优化问题更加复杂。

总而言之，基于非凸变分最小化优化策略的经典光度立体重建方法，通过引入合适的正则化项，有效地解决了传统方法中存在的非凸性、噪声敏感性以及解的不唯一性等问题。然而，如何选择合适的正则化项、设计高效的优化算法、以及处理更复杂的反射模型仍然是该领域的研究热点。未来的研究方向可能包括：开发更鲁棒的算法，处理不同材料和光照条件下的场景，以及将光度立体与其他三维重建技术相结合，从而实现更高精度、更高效率的三维重建。

📣 部分代码

function [func_f,zu,zv,grad_f] = func_f(z,I,rho,M,imask,S,II1,II2,JJ1,JJ2,Ivect,grad_approx)

    

nimgs = size(I,3);

[func_f,A,b,dz_curr,zu,zv] = func_Ab(z(imask),I,rho,M,imask,S,II1,II2,JJ1,JJ2,Ivect);

    % Compute gradient

    if(nargout>=4)

        npix = length(imask);

        

        D = A*M;

        

        if grad_approx

            grad_f = D'*(A*M*z(imask)-b)./(nimgs);

        else

            D_t = D';

            rho_vec = rho(imask);

            M_t = M';

            zgrad = [zu(imask), zv(imask)];

            zgrad_t = zgrad';

            MMz = bsxfun(@times,M_t,(zgrad_t(:))');

            

            MMzSum = MMz(:,1:2:end-1)+MMz(:,2:2:end);

            

            grad = [zgrad, -ones(npix,1)];

            

            tmp1 = (ones(nimgs,1)*(rho_vec./(dz_curr(imask)).^3)').*(S*grad');

            

            tmp2 = cell(npix,1);

            tmp3 = cell(npix,1);

            tmp4 = cell(npix,1);

            for j=1:npix

                tmp2{j,1} = tmp1(:,j)';

            end

            for j=1:npix

                tmp3{j,1} = MMzSum(:,j);

            end

            parfor j=1:npix

                tmp4{j,1} = sparse(tmp3{j,1}*tmp2{j,1});

            end

            D_t = D_t+cell2mat(tmp4');

            

            grad_f = D_t*(A*M*z(imask)-b)./(nimgs);

        end

    end

end

function [func_f,A,b,dz,zu,zv]=func_Ab(z,I,rho,M,imask,S,II1,II2,JJ1,JJ2,Ivect)

    [nrows,ncols,nimgs] = size(I);

    % Depth gradient : 1st order forward diff

    zu = zeros(nrows,ncols);

    zv = zeros(nrows,ncols);

    zu(imask) = M(1:2:end-1,:)*z;

    zv(imask) = M(2:2:end,:)*z;

    % Norm of [-nabla z,1]

    dz2 = zu.^2+zv.^2+1;

    dz = sqrt(dz2);

    % A' field

    rho_over_d = transpose(rho(imask)./dz(imask));

    A1 = -bsxfun(@times,S(:,1),rho_over_d); % mxn

    A2 = -bsxfun(@times,S(:,2),rho_over_d); % mxn

    try

A = sparse2([II1;II2],[JJ1(:);JJ2(:)],[A1(:);A2(:)]);

catch

A = sparse([II1;II2],[JJ1(:);JJ2(:)],[A1(:);A2(:)]);

end

    % b' field

    b = Ivect(imask,:)' - bsxfun(@times,S(:,3),transpose(rho(imask)./dz(imask)));

    b = b(:);

func_f = 0.5*sum((A*M*z-b).^2)/(nimgs);

end

⛳️ 运行结果

🔗 参考文献

🎈 部分理论引用网络文献，若有侵权联系博主删除

👇 关注我领取海量matlab电子书和数学建模资料

🎁  私信完整代码和数据获取及论文数模仿真定制

擅长领域：

🌈 各类智能优化算法改进及应用

生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱调度、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化、CVRP问题、VRPPD问题、多中心VRP问题、多层网络的VRP问题、多中心多车型的VRP问题、 动态VRP问题、双层车辆路径规划（2E-VRP）、充电车辆路径规划（EVRP）、油电混合车辆路径规划、混合流水车间问题、 订单拆分调度问题、 公交车的调度排班优化问题、航班摆渡车辆调度问题、选址路径规划问题、港口调度、港口岸桥调度、停机位分配、机场航班调度、泄漏源定位

🌈 机器学习和深度学习时序、回归、分类、聚类和降维

2.1 bp时序、回归预测和分类

2.2 ENS声神经网络时序、回归预测和分类

2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类

2.4 CNN|TCN|GCN卷积神经网络系列时序、回归预测和分类

2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类

2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类

2.7 ELMAN递归神经网络时序、回归\预测和分类

2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类

2.9 RBF径向基神经网络时序、回归预测和分类

2.10 DBN深度置信网络时序、回归预测和分类

2.11 FNN模糊神经网络时序、回归预测

2.12 RF随机森林时序、回归预测和分类

2.13 BLS宽度学习时序、回归预测和分类

2.14 PNN脉冲神经网络分类

2.15 模糊小波神经网络预测和分类

2.16 时序、回归预测和分类

2.17 时序、回归预测预测和分类

2.18 XGBOOST集成学习时序、回归预测预测和分类

2.19 Transform各类组合时序、回归预测预测和分类

方向涵盖风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、用电量预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断

🌈图像处理方面

图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知

🌈 路径规划方面

旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、 充电车辆路径规划（EVRP）、 双层车辆路径规划（2E-VRP）、 油电混合车辆路径规划、 船舶航迹规划、 全路径规划规划、 仓储巡逻

🌈 无人机应用方面

无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配、无人机安全通信轨迹在线优化、车辆协同无人机路径规划

🌈 通信方面

传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化、水声通信、通信上传下载分配

🌈 信号处理方面

信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化、心电信号、DOA估计、编码译码、变分模态分解、管道泄漏、滤波器、数字信号处理+传输+分析+去噪、数字信号调制、误码率、信号估计、DTMF、信号检测

🌈电力系统方面

微电网优化、无功优化、配电网重构、储能配置、有序充电、MPPT优化、家庭用电

🌈 元胞自动机方面

交通流 人群疏散 病毒扩散 晶体生长 金属腐蚀

🌈 雷达方面

卡尔曼滤波跟踪、航迹关联、航迹融合、SOC估计、阵列优化、NLOS识别

🌈 车间调度

零等待流水车间调度问题NWFSP 、 置换流水车间调度问题PFSP、 混合流水车间调度问题HFSP 、零空闲流水车间调度问题NIFSP、分布式置换流水车间调度问题 DPFSP、阻塞流水车间调度问题BFSP

👇