优化的逐次均值量化变换

连续平均量化变换 (SMQT) 算法是一种非线性变换，它揭示了数据的组织或结构，并去除了增益和偏差等属性。 在一篇名为 The Successive Mean Quantization Transform 的论文中，SMQT 被“应用于语音处理和图像处理”。 当应用于图像时，该算法“可以看作是对图像细节的渐进式关注”。

根据另一篇名为 SMQT-based Tone Mapping Operators for High Dynamic Range Images 的论文，它将产生“具有增强细节的图像”。 该论文描述了将这种变换应用于图像的两种技术。

1. 将 SMQT 应用于每个像素的亮度 - 它描述了从每个像素的 RGB 中获取亮度的公式。

2. 分别在 RGB 的每个通道上应用 SMQT - 分别用于通道 R、G 和 B。

下图显示了两种技术的结果：

通过将第二种技术应用于布鲁因剧院的夜间照片，我们可以看到算法如何逐步放大图像的细节并揭示最初被黑暗隐藏的细节：

在本文中，我们将仔细研究该算法的工作原理，并探索一种简单的优化方法，使该算法在实际图像处理应用中的性能更高。

SMQT 算法

在原始论文中，算法是用抽象的方式描述的。 为了更好地理解 SMQT，我们将介绍一些简单的示例。

假设我们有以下向量（你可以把它想象成一个数组）：

对于彩色图像，我们可以将其视为三个独立的向量，每个向量代表一个特定的颜色通道 (RGB)，向量的每个元素代表相应像素的强度。

应用 SMQT 算法涉及的步骤是：

1. 计算向量的平均值，在本例中为 29.58。

2. 将向量一分为二，将小于或等于 29.58 的数字留在左半部分，将更大的数字留在右半部分：

   15	4	0	5	18	32	48	60	64	59	47	31
   0	0	0	0	0	1	1	1	1	1	1	1

   第二行是我们将为每个值创建的代码。 小于或等于 29.58 的值在其代码中得到 0，大于 29.58 的数字得到 1（这是二进制）。

3. 现在，两个结果向量都按照相同的规则以递归方式单独拆分。 在我们的例子中，第一个向量的平均值是 (15 + 4 + 0 + 5 + 18) / 5 = 8.4，第二个向量的平均值是 (32 + 48 + 60 + 64 + 59 + 47 + 31) / 7 = 48.7。 这两个向量中的每一个都被进一步拆分为另外两个向量，并且根据与平均值的比较，将第二位代码添加到每个向量中。 这是结果：

   4	0	5	15	18	32	47	31	48	60	64	59
   00	00	00	01	01	10	10	10	11	11	11	11

   请注意，对于小于或等于平均值​​的数字（对于每个向量）再次附加 0，对于大于平均值的数字附加 1。

4. 递归地重复该算法：

   0	4	5	15	18	32	31	47	48	60	64	59
   000	001	001	010	011	100	100	101	110	111	111	111

   0	4	5	15	18	31	32	47	48	60	59	64
   0000	0010	0011	0100	0110	1000	1001	1010	1100	1110	1110	1111

   0	4	5	15	18	31	32	47	48	59	60	64
   00000	00100	00110	01000	01100	10000	10010	10100	11000	11100	11101	11110

   此时每个向量只有一个元素。 因此，对于每个连续的步骤，将附加一个 0，因为该数字将始终等于其自身（0 的条件是该数字必须小于或等于均值，即其自身）。

到目前为止，我们的量化级别为 L=5。 如果我们要使用 L=8，我们必须在每个代码后面附加三个 0。 将每个代码从二进制转换为整数（对于 L=8）的结果将是：

最终向量按升序排序。 为了获得输出向量，我们必须将输入向量的原始值替换为其代码。

请注意，在结果向量中：

• 增益被删除。 原来的增益很低，范围从 0 到 64。现在它的增益从 0 到 240，几乎是整个 8 位范围。 据说增益被移除是因为如果我们将原始向量的所有元素乘以一个数字（例如 2），或者如果我们将所有元素除以 2，则输出向量将大致相同。 它的范围大约是目标向量的整个范围（L=8 时为 8 位）。 如果我们将输入向量乘以 2，则输出向量实际上是相同的，因为在每一步中，低于或高于平均值的相同数字将继续低于或高于平均值，我们将添加相同的 0 或 1到输出。 只有当我们除以它时，结果才会大致相同，而不是完全相同，因为像 15 这样的奇数必须四舍五入，然后计算可能会有所不同。 我们使用整个范围从深色图像变为像素范围从深色 (0) 到浅色 (240) 的图像。
• 偏差被消除了，尽管在这个例子中我们不能完全观察到它。 这是因为我们没有偏差，因为我们的最小值是 0。如果我们真的有偏差，它就会被删除。 如果我们取任何数字“K”并将其添加到输入向量的每个元素中，那么每一步中高于和低于平均值的数字的计算将不会发生变化。 输出仍然是相同的。 这也会使太亮的图像无法变成从深色到浅色的图像。
• 我们有一张细节增强的图像。 请注意，对于我们采取的每个递归步骤，我们实际上是如何放大细节，并通过尽可能多地揭示这些细节来构建输出。 由于我们会将这种技术应用于每个 RGB 通道，我们将尽可能多地显示细节，尽管会丢失有关每种颜色的原始色调的信息。

给定一个像 RGB 像素向量的图像，每个点对于 R 来说是 8 位，对于 G 来说是 8 位，对于 B 来说是 8 位； 我们可以从中提取三个向量，每种颜色一个，然后将算法应用于每个向量。 然后我们再次从这三个输出向量中形成 RGB 向量，并且我们拥有经过 SMQT 处理的图像，就像在 Bruin Theatre 中所做的那样。

复杂

该算法要求对于每个量化级别 (L)，必须在第一遍中检查 N 个元素以找到每个向量的均值，然后在第二遍中，必须将这 N 个元素中的每一个与相应的均值进行比较依次拆分每个向量。 最后，N 个元素必须用它们的代码代替。 因此算法的复杂度是O((L*2*N) + N) = O((2*L + 1)*N)，也就是O(L*N)。

从论文中提取的图与 O(L*N) 复杂度一致。 L 越低，以近似线性方式处理的速度越快（每秒的帧数越多）。 为了提高处理速度，论文建议降低 L 的值：“选择较低的 L 级别可以直接降低处理速度，但会降低图像质量”。

改进

在这里，我们将找到一种在不降低图像质量的情况下提高速度的方法。 我们将分析转换过程并找到更快的算法。 为了更全面地了解这一点，让我们看一个重复数字的例子：

向量不能再被分割，并且必须附加零，直到我们得到所需的 L。

为简单起见，假设输入可以从 0 到 31，输出可以从 0 到 7（L=3），如示例中所示。

请注意，计算第一个向量的平均值 (16 + 25 + 31 + 31 + 25 + 16 + 7 + 1 + 1 + 7) / 10 = 16 得出的结果与计算整个最后一个向量的平均值相同，因为它是只是相同的向量，元素的顺序不同：(1 + 1 + 7 + 7 + 16 + 16 + 25 + 25 + 31 + 31) / 10 = 16。

在第二步中，我们必须递归地计算每个向量。 所以我们计算灰色输入的平均值，即前 6 个元素 ((16 + 16 + 7 + 1 + 1 + 7) / 6 = 8)，以及空白输入的平均值，即最后 4 个元素((25 + 31 + 31 + 25) / 4 = 28)：

再次注意，如果我们使用最后一个向量，即完全有序的向量，结果是相同的。 对于前 6 个元素，我们有：((1 + 1 + 7 + 7 + 16 + 16) / 6 = 8)，对于后 4 个元素：((25 + 25 + 31 + 31) / 4 = 28) . 因为它只是一个加法，所以加法的顺序并不重要。

这同样适用于下一步：

计算与最后一个向量相同：((7 + 1 + 1 + 7) / 4 = 4) 将等于 ((1 + 1 + 7 + 7) / 4 = 4)。

这同样适用于其余的向量和步骤。

平均值的计算只是区间中元素的值的总和，除以区间中的元素数。 这意味着我们可以预先计算所有这些值并避免重复此计算 L 次。

让我们看看将我们称之为 FastSMQT 算法的方法应用于我们的示例的步骤：

1. 创建一个包含 32 列和 3 行的表，如下所示。 表中的第一行表示表的索引（0 到 31），在对表进行编码时不必包括在内。 但是它被明确地展示出来以使这个例子更容易理解。

2. 计算每个值的频率后，迭代 N 个元素。 在我们的示例中，元素 1、7、16、25 和 31 的频率均为 2。所有其他元素的频率均为 0。这一步的复杂度为 O(N)。

3. 现在已经填充了频率向量，我们需要迭代该向量（频率向量，而不是输入向量）。 我们不再迭代 N 个元素，而是迭代 R（R 在范围内：2^L），在本例中为 32（0 到 31）。 处理像素时，N 可以是数百万（百万像素），但 R 通常为 256（每种颜色一个向量）。 在我们的示例中，它是 32。我们将同时填充表格的以下两行。 这些行中的第一行（表的第二行）将包含到目前为止的频率总和。 第二个（表的第三个）将包含到目前为止出现的所有元素的值的总和。

   在我们的示例中，当我们到达 1 时，我们将 2 放在表的第二行，因为到目前为止已经出现了两个 1。 我们还在表格的第三行放了一个 2，因为 1 + 1 = 2。我们继续在接下来的每个位置上写入这两个值，直到达到 7。由于数字 7 出现了两次，我们将 2 添加到第二行的累加器，因为到目前为止出现了 4 个数字（1、1、7、7）。 我们将 7 + 7 = 14 添加到第三行，得到 2 + 14 = 16，因为 1 + 1 + 7 + 7 = 16。我们继续这个算法，直到我们完成迭代这些元素。 这一步的复杂度是 O(2^L)，在我们的例子中是 O(32)，当使用 RGB 像素时是 O(256)。

   一世	0	1	2	...	6	7	8	...	15	16	17	...	24	25	26	...	30	31
   n	0	2	0	...	0	2	0	...	0	2	0	...	0	2	0	...	0	2
   n-累积	0	2	2	...	2	4	4	...	4	6	6	...	6	8	8	...	8	10
   和	0	2	2	...	2	16	16	...	16	48	48	...	48	98	98	...	98	160

4. 下一步是从表中删除频率为 0 的元素的列，为了更清楚地查看示例，我们还将删除第二行，因为它不再相关，如您所见。

   一世	1	7	16	25	31
   n	2	4	6	8	10
   和	2	16	48	98	160

5. 请记住，我们可以使用最后一个（完全有序的）向量来为每个步骤进行计算，结果仍然是相同的。 请注意，在此表中，我们有最后一个向量，其中已经进行了所有预先计算。

   我们基本上需要对这个小向量做 SMQT 算法，但与我们开始的原始向量不同，这个会容易得多。

   首先，我们需要计算所有样本的平均值。 我们可以通过 sum[31] / n[31] = 160 / 10 = 16 轻松找到它。所以我们将表格分成两个向量，并开始为每个向量编写二进制代码：

   一世	1	7	16	25	31
   n	2	4	6	8	10
   和	2	16	48	98	160
   代码	0	0	0	1	1

   现在我们再次拆分每个向量。 第一个向量的平均值为：sum[16] / n[16] = 48 / 6 = 8。第二个向量的平均值为：(sum[31] – sum[16]) / (n[31] – n[16]) = (160 – 48) / (10 – 6) = 112 / 4 = 28。

   一世	1	7	16	25	31
   n	2	4	6	8	10
   和	2	16	48	98	160
   代码	00	00	01	10	11

   只剩下一个向量要拆分。 平均值为：sum[7] / n[7] = 16 / 4 = 4。

   一世	1	7	16	25	31
   n	2	4	6	8	10
   和	2	16	48	98	160
   代码	000	001	010	100	110

   如您所见，如果我们遵循原始算法，每个元素的代码都是相同的。 我们在一个比 N 个元素小得多的向量上进行了 SMQT，并且我们还预先计算了我们需要的所有值以找到平均值，因此计算结果向量非常简单快速。

   这一步的复杂度是 O(L*(2^L))，对于 L=8，在实际的图像处理需要中，它是 O(8*(256)) = O(2048) = 常数。

6. 最后一步是迭代 N 个元素，再次用每个元素的代码替换该元素：element[i] = code[i]。 复杂度为 O(N)。 FastSMQT 的复杂度为 O(N) + O(2^L) + O(L*(2^L)) + O(N) = O(2*N) + O((L + 1)*(2 ^L)) = O(N + L*(2^L)) 。

并行性

两种算法都可以并行化，但如果必须转换多个向量，则每个内核运行一个算法会更有效。 例如，在处理图像时，我们可以在不同的内核中处理每个 RGB 通道，而不是并行化这三个计算中的每一个。

为了并行化 SMQT 算法，当我们将一个向量一分为二时，如果有一个内核可用，我们可以在一个新内核中处理每个子向量（实际上一个子向量在当前内核中，另一个在新内核中）。 这可以递归地完成，直到我们达到可用内核的总数。 代码替换值也可以并行进行。

FastSMQT 算法也可以并行化。 第一步，必须将输入向量划分为可用核数 (C)。 必须为每个部分创建一个频率计数向量，并并行填充。 然后我们将这些向量的值添加到一个结果频率向量中。 可以并行化的下一步是最后一步，当我们用它们的代码替换值时。 这可以并行完成。

复杂性比较

原始 SMQT 算法的总复杂度为 O((2*L + 1)*N)，即 O(L*N)。

FastSMQT 的复杂度为 O(N) + O(2^L) + O(L*(2^L)) + O(N) = O(2*N) + O((L + 1)*(2 ^L)) = O(N + L*(2^L))。

给定一个固定的 L，两者的复杂度都变成了 O(N)。 但是对于 FastSMQT，乘以 N 的常数要小得多，并且正如我们将看到的，它在处理时间上有很大的不同。

下图比较了 SMQT 和 FastSMQT 算法在 L=8 时的性能，图像处理就是这种情况。 该图显示了处理 N 个元素所需的时间（以毫秒为单位）。 请注意，对于 L=8，SMQT 的复杂度（使用常数）是 O(17*N)，而对于 FastSMQT，复杂度是 O(2*N + 2304)。

N 越大，SMQT 处理图像所需的时间越长。 另一方面，对于 FastSMQT，增长要慢得多。

对于更大的 N 值，性能差异更加明显：

这里 SMQT 处理某些图像最多需要几秒钟的时间，而 FastSMQT 仍处于“几毫秒”区域内。

即使有多个内核（在本例中为两个），FastSMQT 显然仍优于仅 SMQT：

FastSMQT 不依赖于 L。另一方面，SMQT 线性依赖于 L 的值。例如，当 N = 67108864 时，SMQT 的运行时间随着 L 值的增加而增加，而 FastSMQT 则不：

结论

并非所有优化技术都适用于所有算法，提高算法性能的关键通常隐藏在对算法工作原理的一些了解中。 此 FastSMQT 算法利用了预先计算值的可能性和平均计算的性质。 因此，它允许比原始 SMQT 更快的处理。 速度不受 L 增量的影响，虽然 L 不能比通常的 8 大很多，因为内存使用量是 2^L，对于 8 来说是可以接受的 256（我们频率表中的列数），但是性能增益具有明显的实际效益。

这种速度上的改进使 MiddleMatter 能够在 iOS 应用程序（和 Android 应用程序）中实现该算法，以改进在弱光条件下拍摄的图片和视频。 通过这一改进，可以在应用程序中执行视频处理，否则对于 iOS 设备来说太慢了。

GitHub 上提供了 SMQT 和 FastSMQT 算法的 C++ 代码。