简介针对多光谱成像与压缩编码孔径成像研究方向这份Matlab代码资源提供了从数据立方体获取到图像重建的完整实现框架。资源面向计算机、电子信息工程、数学等专业的大学生可用于课程设计、期末大作业与毕业设计也可作为科研入门参考。压缩包共30个文件包含20个.m主程序与函数脚本辅以.mat数据文件、BMP/PNG图像样本、C语言与mexa64扩展文件以及TXT说明文档结构覆盖成像模拟、数据读取、重建算法与可视化展示等环节。所有代码采用参数化编程关键参数可灵活调整注释明细、思路清晰并附赠案例数据在Matlab 2014/2019a/2024a环境下均可直接运行。已有99人学习下载。通过研读该套代码可掌握压缩编码孔径成像模型构建、数据立方体获取流程、TV范数正则化重建等核心知识并可直接在此基础上扩展算法或复用为实验工具适合需要快速上手或借鉴完整技术路线的学习者。1. 多光谱成像与压缩编码孔径成像为什么需要数据立方体重建从压缩包内容看这套代码的核心不是“拍照”而是“算出来”。多光谱成像中被测场景是一个 x-y-λ 构成的数据立方体探测器永远是二维平面因此单次曝光只能观测到立方体在某方向的积分投影。压缩编码孔径成像CASSI利用随机编码掩模和色散元件把不同波长的空间信息混叠到同一帧测量中再把重建看作解线性反问题。这正好解释了为什么你会看到 Rfuntwist.m 和 RTfuntwist.m 成对出现一个做投影一个做反投影。这份代码适合两类人一类是做压缩感知课程设计或毕业设计的本科生需要一份能跑通再改参数的 Matlab 工程另一类是做多光谱成像的工程师想快速对比 TwIST 与 TV 正则重建效果。包里既有 staticscene1x150.BMP 这样的标准测试图也有 digitalcamerapicture2.jpg 这样的自然图像还有预计算的 mywl.mat 和 Cu10x1x1x1.mat基本可以做到解压后直接运行 RUNME.m 出图。下面从观测模型开始把每个文件在流水线中的角色拆开。2. 从数据立方体到观测模型Rfuntwist 与 RTfuntwist 的矩阵化实现2.1 前向投影的物理过程色散、编码与积分先看观测模型的离散形式。假设数据立方体 f 的尺寸为 n1n2n3其中 n3 是波段数。掩模 M 是 n1n2 的随机二值图案色散使得第 λ 个波长在水平方向平移 s(λ) 个像素。于是探测器上坐标为 (i,j) 的像元强度为 g(i,j) Σ_λ f(i, js(λ), λ) * M(i,j)。写成矩阵乘积 g A fA 就是前向算子。由于 n1n2 和 n3 的乘积很容易超过百万量级显式存储 A 会占用数 GB 内存因此工程实现里永远用函数来代替矩阵乘法。Rfuntwist.m 和 RTfuntwist.m 就是这一对函数。Rfuntwist 接收一个三维立方体输出二维测量RTfuntwist 接收二维测量输出三维立方体。在 TwIST 迭代中每一步需要计算 A f 和 A^T (A f - g)所以这两个函数的执行效率直接决定整个重建的耗时。代码里的 shiftCube.m 负责按波段平移conv2c.m 负责编码掩模与立方体之间的卷积操作它们是 A 的核心组件。2.2 函数边界与转置一致性检查编写这两个函数时最常见的错误是正向与转置的平移方向不一致或者边界处理方式不同导致 A^T 并不是 A 的严格伴随。下面的代码演示了 shiftCube.m 配合累加的前向算子实现% Rfuntwist.m 前向算子示例 function g Rfuntwist(f, mask, shiftList) % f : n1*n2*nL 数据立方体 % mask : n1*n2 编码掩模 % shiftList: 每个波段的色散平移量长度等于 nL [n1, n2, nL] size(f); g zeros(n1, n2); for lambda 1:nL fShifted shiftCube(f(:, :, lambda), shiftList(lambda)); g g fShifted .* mask; end end这里每循环一个波段就把该波段的图像平移 shiftList(lambda) 个像素再用掩模逐像素调制后累加。shiftList 通常由波长向量乘以色散系数得到单位是像素。注意循环顺序不要反过来必须先调制后累加与物理过程一致。对应的转置算子需要把同一个掩模乘到二维测量上再按相反方向平移并放回立方体对应层% RTfuntwist.m 转置算子示例 function f RTfuntwist(g, mask, shiftList) [n1, n2] size(g); nL numel(shiftList); f zeros(n1, n2, nL); for lambda 1:nL f(:, :, lambda) shiftCube(g .* mask, -shiftList(lambda)); end end逻辑说明转置算子的作用是“把二维残差投影回三维空间”所以必须将掩模点乘残差然后反向平移。正向里用了 shift转置里就一定要用 -shift否则迭代不收敛。参数 shiftList 建议用双精度向量而不是整数倍波长索引因为真实色散很难保证整像素对齐如果 shiftCube.m 只支持整型平移可以用 imtranslate 或 interp2 做亚像素插值。为了验证这两个函数是否严格互为转置可以用随机数据做双边检查x randn(32, 32, 8); y randn(32, 32); shiftList 0 : 1 : 7; % 示例 mask double(rand(32, 32) 0.5); lhs sum(sum(sum(RTfuntwist(y, mask, shiftList) .* x))); rhs sum(sum(y .* Rfuntwist(x, mask, shiftList))); fprintf(lhs%.6e, rhs%.6e\n, lhs, rhs);正常时两边应相等误差在 1e-12 量级。如果差异较大首先检查平移方向和掩模是否被重复乘了两次。这个检查在更换掩模或修改边界条件后都应该重跑一遍。提示mask 在生成后应保存下来前向、转置、TwIST 必须共用同一个掩模否则测量与重建不在同一个子空间。2.3 从文件名看测量数据的组织README.txt 之外代码包里还提供了多个数据文件。从命名习惯可以大致推断它们在流水线中的位置文件角色staticscene1x150.BMP高光谱测试场景通常包含多光谱通道堆叠后的灰度表示staticscenex150xdark.BMP暗场景用于估计传感器暗电流或噪声水平digitalcamerapicture2.jpg自然图像可用来生成模拟编码孔径测量Cu10x1x1x1.mat可能是包含色散系数或掩模参数的 10 波段小矩阵mywl.mat预计算的波长下标或光谱响应标定值1.png, 2.png重建结果或测量图的导出文件便于直接对比数字图像文件一般直接由 imread 读入Matlab 会把 uint8 转成 double 后归一化到 [0,1]。这里有个容易踩的坑BMP 和 JPG 的通道数不同如果 datacube2 系列函数期望的是灰度图而读入的是 RGB会在 size(f,3) 处把波段数误判成 3。遇到这种情况要么在读取后调用 imagedatacube2gray.m 先转灰度要么在 RUNME.m 里显式指定 nL。3. 重建算法的核心TwIST 迭代、TV 正则化与 projk 投影3.1 为什么 TwIST 比梯度下降更快CASSI 的逆问题是不适定的A 的列数体素数远大于行数像元数并且采样被噪声污染。直接最小化二范数会产生严重振荡所以目标函数要加正则项min 0.5||Af-g||^2 τ·TV(f)。常见的梯度下降法每次迭代只做一步更新收敛速度慢而 TwISTTwo-step Iterative Shrinkage-Thresholding利用前两次迭代信息做外推收敛速率从 O(1/k) 提升到线性收敛。这和医学超分辨率图像重建里用的快速迭代收缩-阈值思路一致区别在于这里的测量方程换成 CASSI 的 A。如果换成 Matlab 优化工具箱里的 lsqr 或 lsqnonneg 也能算但面对上百万体素时TwIST 的内存占用和收敛速度更有优势。TwISTmod.m 是 TwIST 的局部修改版通常在循环里增加了迭代次数组件和残差早停。TVnormspectralimaging.m 是主入口mycalltoTVnew.m 负责把数据立方体拉直、调用 TwIST、再恢复成三维结构。这一层封装让使用者不需要频繁修改求解器只改测度模型的文件名即可。3.2 TV 正则diffh.m 与 diffv.m全变分正则惩罚相邻像素/体素的差分能保持边缘而平滑噪声。代码里 diffh.m 和 diffv.m 分别计算水平与垂直方向的差分。一个简单的输出如下% TVnormspectralimaging.m 片段 function tvValue TVnormspectralimaging(im, tau) dh diffh(im); % 水平差分大小与 im 相同边界为零 dv diffv(im); % 垂直差分 gradNorm sqrt(dh.^2 dv.^2 eps); tvValue tau * sum(gradNorm(:)); end说明tau 是正则化强度值越大重建图像越光滑但也会抹掉细小光谱特征eps 防止梯度为零时求导无穷大。diffh/diffv 的内部实现如果采用 imfilter 提取边缘边界填充要选 replicate 而不是 circular否则重建图四周会出现上下游通路。光谱维通常没有带宽限制正规的全变分也可以做相邻波段差分。3.3 projk.m 投影约束压缩感知重建经常需要把解限制在某个凸集里最常用的是非负约束和上下界约束。projk.m 做的就是这件事% projk.m 投影到可行域 function x projk(x, lower, upper) x min(x, upper); x max(x, lower); end参数解释lower 和 upper 可以是标量或与 x 同尺寸的矩阵。标量适用于所有体素一致范围矩阵适用于某些区域已知无效的情况。调用时注意顺序先 clamp 上限再 clamp 下限效果等价但遇到 NaN 时会把 NaN 误夹到边界因此进入 projk 之前要确保残差和迭代解没有 NaN。3.4 参数表与调参建议TwIST 求解器需要设置的参数集中在 mycalltoTVnew.m 或 RUNME.m 的头部。一个常见的参数分组如下参数作用建议范围tau / lambda正则项权重控制平滑程度0.005 ~ 0.1先从小值试maxIterTwIST 最大迭代次数200 ~ 2000看收敛曲线tol相对残差变化阈值1e-5 ~ 1e-3mask编码掩模固定随机种子重复实验一致shiftList每个波段的色散量与标定数据匹配单位像素initX初始解零矩阵或 mywl.mat 热启动首次运行时不要直接上 2000 次先用 maxIter200 看重建是否出现明显条带如果条带来自掩模周期性需要打乱掩模像素分布如果条带来自色散偏移过小需要增大相邻波段间距。tau 过大时图像会像油画一样糊tau 过小时噪声会形成颗粒找到中间值通常需要做 35 组对比实验。注意tau 的取值要和数据归一化范围一致。如果测量值被缩放 255 倍tau 也要相应放大否则重建会直接发散。3.5 初始值的坑mywl.mat 的作用mywl.mat 从名字看是预计算的波长列表实际作用可能是给迭代提供一个靠近真实解的启动点。使用好的初始值能将收敛时间缩短一半尤其是当目标和上一次测量的场景相似时。但不要直接把 mywl.mat 当作最终答案它是辅助数据重建真正靠的是 TwIST 和 TV。如果加载 mat 之后发现类与运行平台不符直接换成 zeros(n1,n2,nL) 也能收敛只是速度稍慢。注意 mywl.mat 里的数值可能包含 uint8 类型要先 double 化再参与计算。4. 代码运行与参数设置RUNME、conv2c、shiftCube 实战排错4.1 RUNME.m 的入口逻辑RUNME.m 是这份工程的主入口。通常它会依次完成四件事加载 staticscene1x150.BMP 和 mywl.mat生成或读取掩模调用 Rfuntwist.m 合成二维测量再调用 mycalltoTVnew.m 重建并显示。在终端里运行cd(D:\CASSI_Matlab); clear; close all; RUNME;如果是在 Linux 服务器上用 headless 模式跑需要改成matlab -nodisplay -nosplash -r cd(/home/user/CASSI_Matlab); RUNME; exit参数说明-nodisplay 不启动图形界面-nosplash 跳过启动闪屏最后用 exit 让 Matlab 自动退出。这种方法很适合批量测试不同 tau 值。4.2 矩阵维度不匹配shiftCube 与 conv2c 的调试最常见的运行错误是 Matrix dimensions must agree。原因通常是 f(:, :, lambda) 与 mask 的尺寸不一致或者 shiftCube 在平移后改变了矩阵大小。我的调试顺序是先检查 size(mask) 和 size(imagedatacube2gray(staticscene)) 是否相等再查看 shiftCube 的边界策略。function fshifted shiftCube(f, shift) % 支持亚像素平移的 shiftCube 参考实现 if mod(shift, 1) 0 fshifted circshift(f, [0 shift]); else fshifted imtranslate(f, [0 shift], FillValues, 0); end endcircshift 是循环移位不会缩小矩阵imtranslate 是亚像素插值填充值可以是零或 NaN。转置算子调用 shiftCube 时传 -shift因此 shift 必须是双精度标量。如果这里传入的是负整数circshift 的方向会反过来不会报错但结果错误建议在 Rfuntwist 入口加 assert(shift0 || abs(shift)n2/2) 防止逻辑越界。conv2c.m 的名字表明它可能是二维卷积的循环卷积版本。编码掩模与图像的卷积在频域实现更快fft2(x) .* fft2(mask) 再 ifft2但要注意掩模尺寸需要 padding 到与图像一致否则循环卷积会出现周期性干扰。如果只想用空间卷积就用 conv2 的 same 参数并预先指定边界模式。错误现象可能原因处理方式Error using Matrix dimensions must agreemask 与立方体平面尺寸不一致统一 resize 到相同尺寸Invalid MEX file 或 .mexa64 无法加载平台或 Matlab 版本不同重新编译 ind2rgb8c.c重建结果全是横向条纹色散偏移过小或掩模周期性强增大 shiftList 间距或打乱掩模迭代了很久但残差不动tol 过严或初始化太差改用 mywl.mat 热启动并放宽 tol4.3 mex 文件缺失与编译包内 ind2rgb8c.mexa64 是 C 编译的 mex用于把索引图像快速转换为 RGB。在旧版 Matlab 上可以直接加载但换了 2024a 或 Linux 平台可能报 Invalid MEX file。解决方案是用附带的 C 源文件重新编译mex -setup C mex ind2rgb8c.c编译成功后当前文件夹下会多出 ind2rgb8c.mexa64 或对应的扩展名随后 datacubeplotter 才能调用。如果编译器报缺少头文件可在 mex 命令里加 -I C:\Program Files\MATLAB\R2024a\extern\include。编译产生的二进制文件只在本机 Matlab 版本可用换机器要重新编译。注意在 Windows 上运行 mex -setup 时编译器选择 “Microsoft Visual C” 而不是内置的 LCC否则 Windows 10 以上系统可能出现链接失败。4.4 性能瓶颈与加速技巧循环逐个波段调用 shiftCube 是主要耗时点。一个常用的优化是把所有波段的平移合并成一次矩阵化操作先生成一个 index 矩阵再用 interp2 一次性插值。如果内存足够也可以把 A 显式生成成稀疏矩阵然后直接用最二乘求解但那样就失去了压缩编码孔径使用大立方体时的意义。另一个技巧是在迭代初期用较大的容差加快收敛后期再收紧 tol或者先用小尺寸如 64*64验证算法再扩大到完整分辨率。这与深度学习 Matlab 工具里的“先在子采样数据上调参最后全量训练”思路类似。如果你需要输出中间过程可在 TwISTmod.m 里每 50 次迭代保存一次 x方便画收敛动画或观察掩模伪影的演化。5. 数据立方体的可视化与验证datacubeplotter、ind2rgb8 使用技巧5.1 可视化从立方体到伪彩色重建结果是三维立方体直接看矩阵数值无法判断好坏。包里的 datacubeplotter 是自定义绘图工具通常会配合 imagedatacube2.m 将数据立方体映射为真彩色图像spectrumRGB.m 负责把光谱响应转换为 RGB 显示值。一个典型的调用是cube mycalltoTVnew(measurement, mask); rgb imagedatacube2(cube, mywl); % mywl 对应各波长的中心波长 imshow(rgb); dispCube(cube, montage); % 逐波段拼接显示如果显示器颜色明显偏色可以检查 colorMatchFcn.m 的色度学矩阵是否来自 sRGB 标准。很多压缩重建代码只关心数值最后的色彩映射随便用 parula 替代导致光谱信息无法直观比较。这里应当使用 imagedatacube2gray.m 得到单波段灰度再用 ind2rgb8 搭配制定 colormap 渲染。ind2rgb8c.c 编译出的 mex 比纯 Matlab 实现快一个数量级在 150 个通道的立方体上优势明显。5.2 验证重建质量的方法除了盯图还需要量化验证。假设原始立方体 groundTruth 存在例如从 staticscene1x150.BMP 提取出的 double 矩阵可以计算 PSNR 和 SSIMpeak max(groundTruth(:)); mse mean((groundTruth(:) - cube(:)).^2); psnr 10 * log10(peak^2 / mse); ssimVal ssim(imagedatacube2gray(cube), imagedatacube2gray(groundTruth)); fprintf(PSNR%.2f dB, SSIM%.4f\n, psnr, ssimVal);参数说明peak 取 groundTruth 的动态范围mse 必须把所有体素拉平成向量后再平均ssim 是图像结构相似性需要先转成 2D 灰度图再比较。如果只关心光谱精度可以取某个坐标 (i,j) 的 1*nL 波长曲线与原始曲线画在同一坐标轴里。另一个常用做法是将重建结果经过 Rfuntwist 生成的二次测量与原始测量之间的残差作为数据一致性指标残差应接近噪声水平。如果残差仍有明显结构说明正则项过强或迭代没有收敛。最后的实践技巧用 datacubeplotter 导出 PNG 时先调用 axis off 再 imwrite否则会导入坐标轴边框批量保存时给文件名加上 tau 和迭代次数方便不同参数的结果并排对比。对于 150 波段的立方体不要把每个通道存为单独图片而是合成一个 montage 网格图配合 colorbar 快速检索异常波段。本文还有配套的精品资源点击获取