代码探险家

这个屌丝很懒，什么也没留下！

热门标签

图像处理算法（一） 3D-LUT插值算法之Matlab仿真_3dlut插值

作者：代码探险家 | 2024-07-15 18:01:40

踩

3dlut插值

在图像处理领域中，经常会用到3D-LUT（3D-Look Up Table，三维查找表）来做颜色增强或者色彩转换。3D-LUT是一种处理视频色彩空间转换、非线性色域之间转换的有效手段，简单解释，就是将一组离散数值经过3D-LUT映射后转换为另一组离散数值。其处理效果如下面两个图所示。

图1 原图像

图2 3D-LUT映射后的图像

本文主要就3D-LUT在工程应用领域中的Matlab仿真做简单介绍。

一、3D-LUT概念

3D-LUT（三维查找表）被广泛地在电影工业和显示工业广泛使用，用在不同的显示器、不同的颜色空间之间做色彩校正、显示调节、效果和优化等方面。也就是说，3D-LUT实际上是颜色的转换，是RGB色彩空间的映射，借助3D-LUT的颜色映射，可以对图像/视频的颜色进行调节、校准，借助3D-LUT可以将输入的（R，G，B）像素映射为（R’，G’，B’）。

3D-LUT表中包含了所有可能的（R、G、B）组合，表中的每一组元素包含了R、G、B三个分量的像素值。以8bit位深图像为例，对于一个色彩深度为8bit的RGB图像，3D-LUT表中一共需要有（2^8)^3组元素，每一组元素都包含了3个8bit的值（对应R、G、B三个通道）。在使用这样的一个3D-LUT处理图像时，需要先将这个LUT表存储在系统中，然后根据输入像素流中的每一个像素值，对LUT查表得到相应的元素，作为处理结果输出。这样的一个LUT最小需要占用的内存容量为：（2^8)^3*3*8bit = 402653184 bit = 48M Byte，对软/硬件系统来说就需要用48M Byte大小的内存来存储3D-LUT中的元素。而且，图像的色彩深度（通道位宽）越大，LUT表容量越大，所需要的内存容量也越大。

色彩深度	内存容量
8bit	48M Byte
10bit	3.75G Byte
12bit	288G Byte

例如，当R、G、B颜色的取值范围都是[0，0，0]到[255，255，255]时，如果要将这个映射完全存储起来，我们需要256^3条颜色信息与之对应，这个数据量是极其庞大的，不利于传播与使用。因此，使用采样的方法对LUT进行表示，然后在使用时通过插值处理对图像进行还原，可以有效降低LUT表的容量大小。例如，我们对整个颜色空间进行等分，使用32*32*32个立方体来代表整个RGB色彩空间，LUT中只存储这些等分端点的RGB值，这样就意味着需要有33^3个采样点来表示映射信息，使用33个采样点的意思是，在每个维度上使用33个采样点将256个颜色样本点等分成32个小的立体空间。但是这些点远远不能描述所有的像素值，因此，在使用3D-LUT时，先根据输入像素值找到其位于哪一个小的立方体内，然后再通过插值算法来估计出映射后的像素值。

从下图示意图可以发现，在3D-LUT中，沿着R、G、B方向，像素值是逐渐递增的。

以色彩深度为8bit的3D-LUT表为例，稀疏化策略是：将一个含有（2^8)^3个点的立体空间沿着R、G、B等分成16个小的立方体，每个立方体的顶点视为一个晶格，则可以得到上图所示的尺寸为17*17*17大小的3D-LUT晶格模型。这样一来，3D-LUT表所需要的内存容量就减少为(17^3)*3*8bit ≈ 14.5K Byte，当尺寸为33*33*33时，内存容量为(33^3)*3*8bit ≈ 105.3K Byte，大大降低了对内存资源的依赖。常见的3D-LUT尺寸有17*17*17、33*33*33、65*65*65。

显而易见的是，这样的3D-LUT显然在很多场景下都不适合使用，尤其是像基于FPGA这样的图像处理系统来说，FPGA片内存储器资源有限，往往只有十几K Byte至几M Byte的容量，例如：像ZYNQ Ultrascale+ MPSoC这种比较高端的FPGA芯片，至多也只有36Mbit URAM、35Mbit BRAM、11Mbit Distributed RAM，即使所有的RAM都用来存储LUT表也远远不够。因此，需要一种更轻量的3D-LUT表来解决容量大的问题。一个比较好的方法是对原有的3D-LUT表进行稀疏化，通过等间隔采样的方法，减少3D-LUT表中的数据量。

稀疏化操作有利于降低内存资源消耗，但是，会造成像素映射无法一一对应了。因为它是降采样抽取，输入的RGB值对应的位置往往是某几个点的中间位置，如果简单地取某个近似点的颜色是可行的，但是精度略差一些，通常的做法是计算输入像素值与LUT表中距离最近的晶格点的差值，重新融合计算。因此，在使用稀疏表处理图像时，需要利用3D-LUT表中已知的几个点，通过立体空间插值的方式得到结果，实现达到较高精度近似效果的同时降低系统对内存资源的消耗。比较常见的插值方法有三线性插值、四面体插值、金字塔插值等算法。

插值算法在3D-LUT中扮演着重要的角色，通过对输入像素值和相邻数据点的计算关系进行合理的估计，实现对图像颜色和亮度的调整，同时输出图像更加真实和细腻。不同的插值算法有着不同的特点和适用范围，可以根据具体应用需求选择合适的算法。插值算法的选择和实现对于3D-LUT的效果和性能都有着重要的影响，因此在使用3D-LUT进行图像/视频处理时，插值算法的选择和优化是一个需要重视的问题。

根据插值策略的不同，可以把三维立体空间插值分为三线性插值（八点六面体）、三棱柱插值（六点五面体）、金字塔插值（五点五面体）和四面体插值（四点四面体），这主要是根据切割立方体的不同方式来区分的。

二、3D-LUT的优势

3D-LUT具有一些突出的特性和优势，3D LUT本身是查找表，将本该经过计算得到的映射结果预先存储下来，用于直接读取，避免了对像素点的映射计算，是一种空间换时间的算法，有高响应速度的保证；灵活程度大，带来更好的效果，常规的滤镜是单一映射，通常会将整个图像的所有像素机械地进行映射，例如，一个暖色调的普通滤镜会让图像内原本已经很暖的颜色变的更暖，甚至失真，而3D-LUT则可以实现控制，不同像素值会有不同的非线性的映射效果；支持硬件加速，高响应速度的图像处理通常需要GPU或FPGA进行加速，3D-LUT可以加载为纹理，通过OpenGL进行绘图风格和效果的自定义自由度大、简单LUT文件的制作比较简单，风格和效果可以由富有创意的艺术家们随心所欲地随时创建；用户自定义自由度大，用户可以随时下载或制作符合规范的LUT，在不修改软件的前提下，任意地设置LUT，实现不同的效果；兼容性高，3D-LUT实现的调色可以在任意没有接入调色功能的应用中（如系统应用、三方应用、浏览器、视频流、相机流），无感、无侵入的实现调色，被调色的目标不会出现任何兼容性问题，也没有任何侵入性，跨平台的移植性也很高。

三、插值算法分类

三线性插值（八点六面体）

与双线性插值原理类似，双线性插值是根据二维平面内正方形的4个顶点，进行3次线性插值计算得到新的插值点，而三线性插值是根据三维空间内立方体的8个顶点，进行7次线性插值得到新的插值点。

首先，沿着R轴进行插值得到2个蓝色点的值，利用灰色点（R0,G0,B1）和（R1,G0,B1）得到蓝色点（r,G0,B1）的值，同理可以得到其它3个蓝色点的值；

$V(r,G_{0},B_{0})=V(R_{0},G_{0},B_{0})(1-{\Delta}_{r})+V(R_{1},G_{0},B_{0}){\Delta}_{r}$

$V(r,G_{0},B_{1})=V(R_{0},G_{0},B_{1})(1-{\Delta}_{r})+V(R_{1},G_{0},B_{1}){\Delta}_{r}$

$V(r,G_{1},B_{0})=V(R_{0},G_{1},B_{0})(1-{\Delta}_{r})+V(R_{1},G_{1},B_{0}){\Delta}_{r}$

$V(r,G_{1},B_{1})=V(R_{0},G_{1},B_{1})(1-{\Delta}_{r})+V(R_{1},G_{1},B_{1}){\Delta}_{r}$

然后，沿着G轴进行插值得到2个绿色点的值，利用蓝色点（r,G0,B1）和（r,G1,B1）得到绿色点（r,g,B1）的值，同理可以得到另一个绿色点（r,g,B0）的值；

$V(r,G_{1},B_{1})=V(R_{0},G_{1},B_{1})(1-{\Delta}_{r})+V(R_{1},G_{1},B_{1}){\Delta}_{r}$

$V(r,g,B_{1})=V(r,G_{0},B_{1})(1-{\Delta}_{g})+V(r,G_{1},B_{1}){\Delta}_{g}$

最后，沿着B轴进行插值得到红色点（r,g,b)值，利用两个绿色点（r,g,B1）和（r,g,B0）得到红色点（r,g,b）的值，就是最终的插值结果。

$V(r,g,b)=V(r,b,B_{0})(1-{\Delta}_{b})+V(r,g,B_{1}){\Delta}_{b}$

其中， $\Delta_{r}={r-{R_{0}} \over {R_{1}-R_{0}}} \Delta_{G}={g-{G_{0}} \over {G_{1}-G_{0}}} \Delta_{B}={b-{B_{0}} \over {B_{1}-B_{0}}}$

上述过程，可以用矩阵运算描述：

$V(r,g,b)= VA \Delta = \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}^{T} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & 0 & -1 & 0 & -1 & 0 & 1 & 0\\ 1 & -1 & -1 & 1 & 0 & 0 & 0 & 0\\ 1 & -1 & 0 & 0 & -1 & 1 & 0 & 0\\ -1 & 1 & 1 & -1 & 1 & -1 & -1 & 1 \end{Bmatrix}^{T} \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \Delta_{b}\Delta_{r}\\ \Delta_{r}\Delta_{g}\\ \Delta_{g}\Delta_{b}\\ \Delta_{r}\Delta_{g}\Delta_{b}\\ \end{Bmatrix}$

实际上，VA的值并不依赖于（r,g,b）值，所以在实际应用时，可以事先把每一个子立方体对应的VA值计算出来存储在内存中，根据输入的（r,g,b）值进行索引，有利于降低插值过程中的计算开销。

不可避免的是，在对每一个像素插值计算时，都需要对3D-LUT表查找8次。对于FPA实现来说，可以采用串行查表或并行查表，这两种策略都有各自的优缺点：

串行查表：使用单个RAM存储3D-LUT表，每个像素需要用8个周期完成查表操作，插值处理难以流水线化，处理速度极慢；但有利于节省RAM存储器资源。
并行查表：使用8个RAM存储3D-LUT表，每个像素可以在一个周期完成查表操作，插值处理容易流水线化，处理速度较快；但RAM存储器资源会提升到8倍。

显然，无论是串行查表策略还是并行查表策略都不是很适合在FPGA这样资源有限的硬件上使用。因此，加快计算速度的一个有效策略是使用立方体的一部分顶点进行插值处理，这是可以达到加快速度、节省RAM资源的一种折中方法。只有三种方法可以将一个立方体分割成多个顶点数量相同的3D结构：三棱柱（6个顶点）、金字塔（5个顶点）、四面体（4个顶点）。

三棱柱插值（六点五面体）

三棱柱插值的思想是将八点六面立方体立方体分割为2个三棱柱，然后判断输入像素值位于其中的哪一个三棱柱中，利用该三棱柱的2个底面各插值出一个点，然后再利用这两个点做线性插值得到新的插值点。

将八点六面立方体分割为2个三棱柱（p1、p2），根据约束条件确定输入像素（r,g,b）位于哪一个三棱柱，约束条件：若 $\Delta b > \Delta r$ ，则选择p1，否则选择 p2；
利用三棱柱的上、下2个底面，做三角插值得到两个插值点；

${V(r,g,b)}^0 = V(R_{0},G_{0},B_{0}){{S_{000}} \over {S}}+V(R_{0},G_{0},B_{1}){{S_{001}} \over {S}}+V(R_{0},G_{0},B_{0}){{S_{101}} \over {S}}$

${V(r,g,b)}^1 = V(R_{0},G_{1},B_{0}){{S_{010}} \over {S}}+V(R_{0},G_{1},B_{1}){{S_{011}} \over {S}}+V(R_{1},G_{1},B_{1}){{S_{111}} \over {S}}$

S_{000} = S_{010}，S_{001} = S_{011}， S_{101} = S_{111}

利用得到的上下底面的2个点做线性插值得到插值点的像素值；

$V(r,g,b)={V(r,g,b)}^0(1-\Delta_{g})+{V(r,g,b)}^1\Delta_{g}$

上述过程，可以用矩阵运算描述：

$V(r,g,b)_{p1}= \Delta B_{1}V = \begin{Bmatrix} 1 & \Delta_{b} & \Delta_{r} & \Delta_{g}&\Delta_{bg} & \Delta_{rg} \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & -1 & 0 & 1 & 0\\ -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & -1 & -1 & 1 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

$V(r,g,b)_{p2}= \Delta B_{2}V = \begin{Bmatrix} 1 & \Delta_{b} & \Delta_{r} & \Delta_{g}&\Delta_{bg} & \Delta_{rg} \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0\\ -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & -1 & 0 & 0 & -1 & 1\\ 1 & -1 & -1 & 1 & 0 & 0 & 0 & 0 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

从上述计算公式可以发现， $B_{1}V$ 和 $B_{2}V$ 的值不依赖于（r,g,b）值，所以在实际应用时，可以事先计算出来存储在内存中，然后根据输入的（r,g,b）值进行索引，这样可以加快插值处理的速度、节省计算资源。

金字塔插值（五点五面体）

金字塔插值的思想是将八点六面立方体分割成3个金字塔（五点五面体），然后判断输入像素值（r,g,b）位于哪一个金字塔中，进而利用这个金字塔插值得到新的插值点。

将八点六面立方体分割为3个金字塔（p1、p2、p3），根据约束条件确定输入像素（r,g,b）位于哪一个金字塔，约束条件：

$\Delta_{r} > \Delta_{b}$ ， $\Delta_{g} > \Delta_{b}$ ，则选择p1;

$\Delta_{b} > \Delta_{r}$ ， $\Delta_{g} > \Delta_{r}$ ，则选择p2。

假设输入像素（r,g,b）位于p1，则可以得到插值结果：

${V(r,g,b)} = V(R_{0},G_{0},B_{0}){{V_{000}} \over {V_{p1}}} +V(R_{0},G_{0},B_{1}){{V_{001}} \over {V_{p1}}} +V(R_{1},G_{0},B_{1}){{V_{101}} \over {V_{p1}}}+V(R_{0},G_{1},B_{1}){{V_{011}} \over {V_{p1}}} +V(R_{1},G_{1},B_{1}){{V_{111}} \over {V_{p1}}}$

$V_{p1}=V_{000}+V_{001}+V_{101}+V_{011}+V_{111}$

上述计算可以用矩阵运算描述：

$V(r,g,b)_{p1}= \Delta_{p1}C_{1}V = \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \Delta_{r}\Delta_{g} \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1\\ -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & -1 & -1 & 1 & 0 & 0 & 0 & 0 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

$V(r,g,b)_{p2}= \Delta_{p2}C_{2}V = \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \Delta_{g}\Delta_{b} \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & 1\\ -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & -1 & 0 & 0 & -1 & 1 & 0 & 0 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

$V(r,g,b)_{p3}= \Delta_{p3}C_{3}V = \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \Delta_{b}\Delta_{r} \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1\\ 1 & 0 & -1 & 0 & -1 & 0 & 1 & 0 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

同样地， $C_{1}V$ 、 $C_{2}V$ 和 $C_{3}V$ 的值不依赖于（r,g,b）值，所以在实际应用时，可以事先计算出来存储在内存中，然后根据输入的（r,g,b）值进行索引，这样可以加快插值处理的速度、节省计算资源。

四面体插值（四点四面体）

四面体插值（Tetrahedral interpolation）的思想是将八点六面立方体分割成6个四点四面体，然后判断输入像素值（r,g,b）位于哪一个四面体中，进而利用这个四面体的4个顶点插值得到新的插值点，而且，用于插值的4个顶点距离插值点的距离最近，所以精度比其他插值算法也更好。

将八点六面立方体分割成6个四面体（t1、t2、t3、t4、t5、t6），根据约束条件确定输入像素（r,g,b）位于哪一个四面体，约束条件：

$\Delta_{b} > \Delta_{r} > \Delta_{g}$ ，则选择t1;

$\Delta_{b} > \Delta_{g} > \Delta_{r}$ ，则选择t2;

$\Delta_{g} > \Delta_{b} > \Delta_{r}$ ，则选择t3;

$\Delta_{r} > \Delta_{b} > \Delta_{g}$ ，则选择t4;

$\Delta_{r} > \Delta_{g} > \Delta_{b}$ ，则选择t5;

否则，选择t6。

得到合适的四面体后，即可按照公式利用对应四面体的4个定点进行插值计算；

$V_{a^{'}}=V(R_{0},G_{0},B_{0})(1-\Delta_{g})+V(R_{1},G_{1},B_{1})\Delta_{g}$

$V_{b^{'}}=V(R_{0},G_{0},B_{1})(1-\Delta_{g})+V(R_{1},G_{1},B_{1})\Delta_{g}$

$V_{c^{'}}=V(R_{1},G_{0},B_{1})(1-\Delta_{g})+V(R_{1},G_{1},B_{1})\Delta_{g}$

$V_{a^{''}}=V{a^{'}}(1-\Delta_{r})+V_{c^{'}}\Delta_{r}$

$V_{b^{''}}=V{b^{'}}(1-\Delta_{r})+V_{c^{'}}\Delta_{r}$

$V_{x}=V{a^{''}}(1-\Delta_{b})+V_{b^{''}}\Delta_{b}$

对于四面体插值，计算过程可以用矩阵运算描述：

$V(r,g,b)_{t1}= \Delta_{t}T_{1}V = \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & -1 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

$V(r,g,b)_{t2}= \Delta_{t}T_{2}V = \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & 1\\ 0 & 0 & 0 & 0 & -1 & 1 & 0 & 0 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

$V(r,g,b)_{t3}= \Delta_{t}T_{3}V = \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & 1\\ -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

$V(r,g,b)_{t4}= \Delta_{t}T_{4}V = \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0\\ -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

$V(r,g,b)_{t5}= \Delta_{t}T_{5}V = \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1\\ -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

$V(r,g,b)_{t6}= \Delta_{t}T_{6}V = \begin{Bmatrix} 1\\ \Delta_{b}\\ \Delta_{r}\\ \Delta_{g}\\ \end{Bmatrix} \begin{Bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1\\ 0 & -1 & 0 & 1 & 0 & 0 & 0 & 0\\ -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \end{Bmatrix} \begin{Bmatrix} V(R_{0},G_{0},B_{0})\\ V(R_{0},G_{1},B_{0})\\ V(R_{1},G_{0},B_{0})\\ V(R_{1},G_{1},B_{0})\\ V(R_{0},G_{0},B_{1})\\ V(R_{0},G_{1},B_{1})\\ V(R_{1},G_{0},B_{1})\\ V(R_{1},G_{1},B_{1}) \end{Bmatrix}$

同样地， $T_{i}V$ 的值不依赖于（r,g,b）值，但是，如果事先将 $T_{i}V$ 的值事先计算好保存在存储器中，意味着需要存储6个四面体对应的 $T_{i}V$ 表，对于在PC端应用来说是没有问题的，确实能起到加速计算的效果，但是对于FPGA这种存储资源有限的器件来说，显然不可取。

四、Matlab仿真

以下仿真是从FPGA实现本算法的角度出发，所有计算过程均作了量化处理，并对插值过程进行了步骤分解，以方便快速实现FPGA移植。


%main.m
 
clear all; 
clc;
%{
工程注释：3D-LUT插值算法性能验证，主要实现了四面体插值与三线性插值两种算法。
1、四面体插值思想：首先，需要根据输入像素的R、G、B三个通道的像素值计算插值权重（delta_r、delta_g 、delta_b），以及3D-LUT表中最靠近插值点的像素坐标BGR_000（index_r、index_g、index_b）；
                   然后，根据3个通道的插值权重大小顺序，选择对应的四面体，根据BGR_000的坐标，在3D-LUT表中查找到该四面体的4个顶点对应的RGB像素值；
                   最后，根据相应的四面体插值公式，利用4个顶点的RGB像素值计算插值结果。
四面体插值公式：
四面体              判断条件                                                   插值公式
  T0    delta_g >= delta_b >= delta_r       (1 - delta_g) * BGR_000 + (delta_g - delta_b) * BGR_010 + (delta_b - delta_r) * BGR_110 + delta_r * BGR_111 
  T1    delta_b >  delta_r >  delta_g       (1 - delta_b) * BGR_000 + (delta_b - delta_r) * BGR_100 + (delta_r - delta_g) * BGR_101 + delta_g * BGR_111 
  T2    delta_b >  delta_g >= delta_r       (1 - delta_b) * BGR_000 + (delta_b - delta_g) * BGR_100 + (delta_g - delta_r) * BGR_101 + delta_r * BGR_111 
  T3    delta_r >= delta_g >  delta_b       (1 - delta_r) * BGR_000 + (delta_r - delta_g) * BGR_001 + (delta_g - delta_b) * BGR_011 + delta_b * BGR_111 
  T4    delta_g >  delta_r >= delta_b       (1 - delta_g) * BGR_000 + (delta_g - delta_r) * BGR_010 + (delta_r - delta_b) * BGR_011 + delta_b * BGR_111 
  T5    delta_r >= delta_b >= delta_g       (1 - delta_r) * BGR_000 + (delta_r - delta_b) * BGR_001 + (delta_b - delta_g) * BGR_101 + delta_g * BGR_111 
2、三线性插值思想：
三线性插值公式：首先，需要根据输入像素的R、G、B三个通道的像素值计算插值权重（delta_r、delta_g 、delta_b），以及3D-LUT表中最靠近插值点的像素坐标 BGR_000（index_r、index_g、index_b）；
                然后，利用BGR_000的坐标，利用BGR_000的坐标计算以BGR_000为原点的正方体的8个顶点坐标，并在3D-LUT表中查找到这8个顶点对应的RGB像素值；
                最后，利用这8个像素值分别按照R、G、B三个方向进行三次线性插值得到插值结果。
三线性插值公式：
    
P_00 = BGR_000 * (1 - delta_r) + BGR_001 * delta_r
P_01 = BGR_100 * (1 - delta_r) + BGR_101 * delta_r
P_10 = BGR_010 * (1 - delta_r) + BGR_011 * delta_r
P_11 = BGR_110 * (1 - delta_r) + BGR_111 * delta_r
P_0 = P_00 * (1 - delta_g) + P_01 * delta_g
P_1 = P_10 * (1 - delta_g) + P_11 * delta_g
 
P = P_0 * (1 - delta_b) + P_1 * delta_b
注：delta_r、delta_g、delta_b表示插值权重，值域为[0,1]；但是，这里是按照FPGA实现思路做的，因此将所有的权重量化成了WEIGHT_QN bit的整数。
3、结论：
    三线性插值与四面体插值实现的效果都非常接近，使用线性表时，两者几乎无差别，有少量的像素差值为1；使用非线性表时，两者之间的差别略大，但也只差值为1的像素数量变多了一些；
%}
%像素通道位宽
PIXEL_CHANNEL_WIDTH = 8;
 
%插值权重量化因子
WEIGHT_QN = 12;
W_Q = 2 ^ WEIGHT_QN - 1;
 
%LUT表量化位宽
LUT_QN = 8;
LUT_Q = 2 ^ LUT_QN - 1;
 
%LUT维度
LUT_DIM = 33;
%{
    第一步：预处理
            对3D-LUT数据做整型量化操作
            对图像数据做预处理
%}
%加载3D-LUT参数
 
lut = load('input_lut_33.txt','%f10');  %非线性表
 
% lut = load('IdentityLUT33.txt','%f');  %线性表
 
%将float型lut参数量化为8bit整数
lut_q = round(lut * LUT_Q);
 
%加载图像数据
input_image = imread('0.bmp');
[row,col,chan]=size(input_image);
%{
    第二步：插值处理
%}
tetrahed_out = [0 0 0];
trilinear_out = [0 0 0];
tetrahed_image = zeros(row,col,chan);
trilinear_image = zeros(row,col,chan);
 
%1、计算LUT坐标和权重
 
for i = 1 : 1: row
    for j = 1 : 1 : col 
        %计算float型坐标
        float_index_r = double(input_image(i,j,1)) * double(LUT_DIM - 1) / double(255);
        float_index_g = double(input_image(i,j,2)) * double(LUT_DIM - 1) / double(255);
        float_index_b = double(input_image(i,j,3)) * double(LUT_DIM - 1) / double(255);        
           
        %得到整型坐标
        index_r = floor(float_index_r);
        index_g = floor(float_index_g);
        index_b = floor(float_index_b);
         
        %边界条件判断 边界约束     
 
        if index_r == LUT_DIM - 1
            index_r = index_r - 1; 
        end
        if index_g == LUT_DIM - 1
            index_g = index_g - 1;
        end
        if index_b == LUT_DIM - 1
            index_b = index_b - 1;
        end
        
        
        %计算权重 量化 放大 2^Q_N 倍
        weight_r = floor((double(float_index_r) - double(index_r)) * double(W_Q));
        weight_g = floor((double(float_index_g) - double(index_g)) * double(W_Q));
        weight_b = floor((double(float_index_b) - double(index_b)) * double(W_Q));
        
 
        %计算LUT中8个像素的地址
        
        index_000 = index_b * LUT_DIM * LUT_DIM + index_g * LUT_DIM + index_r + 1;                    % BGR 000
        index_001 = index_b * LUT_DIM * LUT_DIM + index_g * LUT_DIM + index_r + 1 + 1;                % BGR 001
        index_010 = index_b * LUT_DIM * LUT_DIM + (index_g + 1) * LUT_DIM + index_r + 1;              % BGR 010
        index_011 = index_b * LUT_DIM * LUT_DIM + (index_g + 1) * LUT_DIM + index_r + 1 + 1;          % BGR 011
        index_100 = (index_b + 1) * LUT_DIM * LUT_DIM + index_g * LUT_DIM +  index_r + 1;             % BGR 100
        index_101 = (index_b + 1) * LUT_DIM * LUT_DIM + index_g * LUT_DIM + index_r + 1 + 1;          % BGR 101
        index_110 = (index_b + 1) * LUT_DIM * LUT_DIM + (index_g + 1) * LUT_DIM + index_r + 1;        % BGR 110
        index_111 = (index_b + 1) * LUT_DIM * LUT_DIM + (index_g + 1) * LUT_DIM + index_r + 1 + 1;    % BGR 111      
                    
        %LUT寻址
        
        if weight_g >= weight_b && weight_b >= weight_r     %t1
            %LUT表中4个坐标的下标
            for k = 1 : 1 : 3
                tetrahed_out(k) = tetrahed_interpl(lut_q(index_000,k),lut_q(index_010,k),lut_q(index_110,k),lut_q(index_111,k),weight_g,weight_b,weight_r,W_Q);
            end
       
        elseif weight_b > weight_r && weight_r > weight_g   %t2
        %LUT表中4个坐标的下标  
            for k = 1 : 1 : 3
                tetrahed_out(k) = tetrahed_interpl(lut_q(index_000,k),lut_q(index_100,k),lut_q(index_101,k),lut_q(index_111,k),weight_b,weight_r,weight_g,W_Q);
            end
            
        elseif weight_b > weight_g && weight_g >= weight_r   %t3
            %LUT表中4个坐标的下标
            for k = 1 : 1 : 3
                tetrahed_out(k) = tetrahed_interpl(lut_q(index_000,k),lut_q(index_100,k),lut_q(index_110,k),lut_q(index_111,k),weight_b,weight_g,weight_r,W_Q);
            end
            
        elseif weight_r >= weight_g && weight_g > weight_b   %t4
            %LUT表中4个坐标的下标         
            for k = 1 : 1 : 3
                tetrahed_out(k) = tetrahed_interpl(lut_q(index_000,k),lut_q(index_001,k),lut_q(index_011,k),lut_q(index_111,k),weight_r,weight_g,weight_b,W_Q);
            end
            
        elseif weight_g > weight_r && weight_r >= weight_b   %t5   
            %LUT表中4个坐标的下标
            for k = 1 :  3 %R G B三个通道插值
                tetrahed_out(k) = tetrahed_interpl(lut_q(index_000,k),lut_q(index_010,k),lut_q(index_011,k),lut_q(index_111,k),weight_g,weight_r,weight_b,W_Q);
            end 
                      
         elseif weight_r >= weight_b && weight_b >= weight_g   %t6
            %LUT表中4个坐标的下标
            for k = 1 : 1 : 3
                tetrahed_out(k) = tetrahed_interpl(lut_q(index_000,k),lut_q(index_001,k),lut_q(index_101,k),lut_q(index_111,k),weight_r,weight_b,weight_g,W_Q);
            end
            
        end
 
       tetrahed_image(i,j,:) = tetrahed_out;
        
        %三线性插值
        for k = 1 : 1 : 3
            trilinear_out(k) = trilinear_interpolation(lut_q(index_000,k),lut_q(index_001,k),lut_q(index_010,k),lut_q(index_011,k),lut_q(index_100,k),lut_q(index_101,k),lut_q(index_110,k),lut_q(index_111,k),(W_Q - weight_r),(W_Q - weight_g),(W_Q - weight_b),W_Q);          
        end
        trilinear_image(i,j,:) = trilinear_out;
    end
end
%保存输出图像
 
tetrahed_image = tetrahed_image ./ (2 ^ (LUT_QN - PIXEL_CHANNEL_WIDTH + WEIGHT_QN));
trilinear_image = trilinear_image ./ (2 ^ (LUT_QN - PIXEL_CHANNEL_WIDTH + WEIGHT_QN * 3));
 
tetrahed_image = floor(tetrahed_image);
trilinear_image = floor(trilinear_image);
 
tetrahed_image(1,1,:);
 
imwrite(uint8(tetrahed_image),'tetrahed_image.bmp');
imwrite(uint8(trilinear_image),'trilinear_image.bmp');
 
 
figure(1);
%显示输入图像
subplot(2,2,1);
imshow(uint8(input_image));
title('原始输入图像');
 
%显示四面体插值图像
subplot(2,2,2);
imshow(uint8(tetrahed_image));
title('Matlab实现四面体插值图像');
 
%显示三线性插值图像
subplot(2,2,3);
imshow(uint8(trilinear_image));
title('Matlab实现三线性插值图像');
 
 
disp('插值权重量化位宽');
disp(WEIGHT_QN);
 
disp('LUT表量化位宽');
disp(LUT_QN);
 
diff = abs(tetrahed_image - trilinear_image);
 
R_Error_Sum = 0;
R_Error_Num = 0;
G_Error_Sum = 0;
G_Error_Num = 0;
B_Error_Sum = 0;
B_Error_Num = 0;
for i = 1 : 1 : row
    for j = 1 : 1 : col
        if(diff(i,j,1) > 1)
            disp('Error!');
            disp(diff(i,j,1));
            R_Error_Num = R_Error_Num + 1;
            R_Error_Sum = R_Error_Sum + diff(i,j,1);
        end 
        if(diff(i,j,2) > 0)
            disp('Error!');
            disp(diff(i,j,2));
            G_Error_Num = G_Error_Num + 1;
            G_Error_Sum = G_Error_Sum + diff(i,j,2);
        end 
        if(diff(i,j,3) > 0)
            disp('Error!');
            disp(diff(i,j,3));
            B_Error_Num = B_Error_Num + 1;
            B_Error_Sum = B_Error_Sum + diff(i,j,3);
        end 
    end 
end 
disp('R_Error_Sum is');
disp(R_Error_Sum);
 
disp('R_Error_Num is');
disp(R_Error_Num);
 
disp('G_Error_Sum is');
disp(G_Error_Sum);
 
disp('G_Error_Num is');
disp(G_Error_Num);
 
disp('B_Error_Sum is');
disp(B_Error_Sum);
 
disp('B_Error_Num is');
disp(B_Error_Num);
 
Avg_r = R_Error_Sum / R_Error_Num
Avg_g = G_Error_Sum / G_Error_Num
Avg_b = B_Error_Sum / B_Error_Num
 
 
%%%%%%%%%%%%        Matlab实现的三线性插值图像与Matlab实现的四面体插值图像对比        %%%%%%%%%%%%
 
MES = sum(sum((double(tetrahed_image) - double(trilinear_image)).^2))/(row * col);     
PSNR = 20 * log10(255 / sqrt(MES));  
 
disp('Matlab实现的三线性插值图像与Matlab实现的四面体插值图像性能对比');
 
disp('R通道的峰值信噪比：');
disp(PSNR(1));
disp('G通道的峰值信噪比：');
disp(PSNR(2));
disp('B通道的峰值信噪比：');
disp(PSNR(3));
 
disp('R通道的均方误差：');
disp(MES(1));
disp('G通道的均方误差：');
disp(MES(2));
disp('B通道的均方误差：');
disp(MES(3));
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
 
 
%四点四面体差值
function pixel = tetrahed_interpl(p0,p1,p2,p3,w0,w1,w2,len)
    
    pixel = (len - w0) * p0 + (w0 - w1) * p1 + (w1 - w2) * p2 + w2 * p3;
    
end
 
%三线性插值
function pixel = trilinear_interpolation(p1,p2,p3,p4,p5,p6,p7,p8,delta_r,delta_g,delta_b,len)
    p_00 = linear_intpl(p1,p2,delta_r,len);
    p_01 = linear_intpl(p3,p4,delta_r,len);
    p_10 = linear_intpl(p5,p6,delta_r,len);
    p_11 = linear_intpl(p7,p8,delta_r,len);
        
    p_0 = linear_intpl(p_00,p_01,delta_g,len);
    p_1 = linear_intpl(p_10,p_11,delta_g,len);
        
    pixel = linear_intpl(p_0,p_1,delta_b,len);
end

结论：从峰值信噪比（PSNR）指标来看，四面体插值得到的图像和三线性插值得到的图像，各通道的PSNR都在50db以上，但是由于四面体插值时做了量化处理，因此和三线性插值存在一定的误差，各通道像素误差在2以内，基本上达到了想要的效果。

附：Color Space Conversion

总结

本文主要是对3D-LUT插值算法原理做了简单的介绍，并且利用Matlab对算法进行了仿真，为后续FPGA移植打下了坚实的理论基础。

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/代码探险家/article/detail/830500