因为论文太贵了，所以我照着论文标题和摘要尝试复现了一篇论文。
上一篇文章讨论了DFT隐写，这种隐写方法能够抵抗各种攻击，但总体来说那更像是盲水印，用于标记而非传递信息。这一篇文章本来想写DCT，但后来我发现单纯的DWT和DCT隐写都是在频域上做手脚，并且效果都不是特别好。后来我看到一篇论文讨论DWT-DCT-SVD隐写，但很遗憾，这篇论文收费的，我在sci-hub上也没有找到，所以就照着论文标题和摘要，问了问ChatGPT，然后基于我和他的解读写了这篇文章。
这篇文章并没有前几篇写的那么详尽，因为我对这之中提及的概念并不十分熟练，只是现学现卖。同时这种隐写方法的表现也没有那么好，或者说需要大量的尝试才能发挥出最佳水平。因此这篇文章描述的流程可能和论文里完全不一样，也可能是错的，也可能我漏掉了一些细节，总之我不太想在这种不确定（不清晰）的路径上花费太多精力去打磨。所以这篇文章可能在质量上不如前两篇。
论隐写容量还得是LSB啊
参考论文
这里参考的论文是DWT-DCT-SVD based watermarking。摘要如下：
Some works are reported in the frequency domain watermarking using Single Value Decomposition (SVD). The two most commonly used methods are based on DCT-SVD and DWT-SVD. The commonly present disadvantages in traditional watermarking techniques such as inability to withstand attacks are absent in SVD based algorithms. They offer a robust method of watermarking with minimum or no distortion. DCT based watermarking techniques offer compression while DWT based compression offer scalability. Thus all the three desirable properties can be utilized to create a new robust watermarking technique. In this paper, we propose a method of non-blind transform domain watermarking based on DWT-DCT-SVD. The DCT coefficients of the DWT coefficients are used to embed the watermarking information. This method of watermarking is found to be robust and the visual watermark is recoverable without only reasonable amount of distortion even in the case of attacks. Thus the method can be used to embed copyright information in the form of a visual watermark or simple text.
隐写方案
我拿着论文摘要去问ChatGPT，可能是因为我没开会员，所以ChatGPT也没说出来个所以然。虽然是2008年的论文，可能因为没有公开访问，所以ChatGPT的训练资料中没有相关的内容，总之帮助不大。从论文摘要来看，首先需要将图片应用DWT，然后对DWT结果再进行DCT，然后对DCT结果应用SVD，最后在SVD的S矩阵里进行隐写。
DWT是离散小波变换；DCT是离散余弦变换，SVD是线性代数中的奇异值分解。
话是这么说，在写DFT隐写的时候我就注意到了，除非有诸如FFT这种经过优化的快速算法，否则按照原始的公式定义去处理大图片会非常慢。但是就我在网络上搜索的结果来看，好像DWT和DCT都没有什么特别好用的快速算法，并且我期待这种隐写结果能够抵御JPEG压缩，所以我认为最好的办法就是从JPEG压缩中找灵感。
JPEG压缩
JPEG压缩是利用人眼对高频信息和色彩不敏感的原理，丢弃图像中不易被人眼辨识的信息，从而达到减少文件大小的效果。总地来说，JPEG压缩的流程如下：
- 将图像从RGB空间转换成YUV空间，其中Y（明度）完全保留，U和V进行下采样（即2x2的格子只取左上角的值，相当于横竖各砍掉一半）
- 以8x8像素划分成区块（Y通道是8x8，U和V对应是4x4），对每个区块计算DCT。DCT能够将能量聚集到矩阵的上三角部分（左上是低频，右下是高频）
- 对DCT运算结果进行量化（取整）。根据量化表对不同频率的值进行整除。例如低频量化值是2，那么相对保留的细节就多一些（可取的值就是2的倍数），高频量化值相对高，例如99，则细节就少（可取的值只有99的倍数）。经过量化后能够做到下三角部分（高频）几乎都为0
- 对量化结果进行编码
这个过程中的主要信息损失就在量化部分。DCT也有损失，但只是损失了小数点，而量化则是以倍数进行损失的。例如量化值是2，那么5就会被量化成4；如果量化值是99，那么0到98都会被量化成0，而99到197全都会被量化成1，这种损失是非常大的。此外JPEG压缩是以8x8像素为单位独立压缩的，所以我们也可以把图片分成8x8的像素块进行处理。
DWT
离散小波变换是小波变换的离散形式（这不是废话吗）。什么是小波变换呢？本质上它是一组滤波器。在分解的时候，把一个信号输入到低通滤波器和高通滤波器，分别得到这个信号的低频分量和高频分量。但是这种分量并不是完全分离的，例如在f/2附近，低频分量可能会包含一些这个频率，而高频分量也会包含一些这个频率。对于这两个分量来说，因为频率被减半了，因此根据采样定理，这些分量可以被下采样（subsampled by 2）而不会丢失信息。在重组的时候，我们将对应的高低频分量输入到对应的滤波器中，然后通过插值还原出原始信号。
至于为什么采样定理说可以通过下采样丢弃样本而不丢失信息，以及为什么通过插值就能还原出原始信号，我真没学过信号处理，我真的不知道。就上面这点还是我前几周刚学的呢。
但是眼尖的读者可能会发现一些问题：你刚才说的什么信号啊滤波器啊，听起来都是一维的，图片是二维，那怎么办呢？说实话，我也在这个问题上反复了好久，看了好多公式和代码，就是没能理解为什么要那样处理。最后还是维基百科一个图片解决了我的疑惑。因为我懒得下载图片了，这里口述一下：
图片作为原始信号，如果我们以行为单位，那么每行就是一个一维信号。这时将这个一维信号进行DWT，得到的低频是原始行信号的近似，而高频则是水平方向的突变（竖直边）。我们将这个低频信号记作L，高频信号记作H。这时L和H是等高半宽，因为按行（沿水平方向）滤波使得宽度减半，而高度没有变。我们可以很自然地将L放在左边，H放在右边，这样我们就得到了和原来大小一样的矩阵了。
这个时候再按列看，对于L来说，将每一列作为一维信号进行DWT，得到了LL（低频）和LH（高频）。其中LL经过了行和列两个方向的低通滤波，因此留下的是原始图片的低频信息，这个矩阵也被称作近似矩阵。而LH经过按列的高通滤波，留下的是原始图像的竖直方向的突变，也就是水平边。
类似地，如果对H也进行按列的DWT，会得到HL和HH。其中HL是对竖直边H的近似，因此留下的还是竖直边的特征。对于HH由于是对竖直边H进行竖直方向的滤波，留下的是水平方向的边，而原始信号H已经是竖直边了，综合来看HH矩阵保留的是原始图片的斜边特征。
按列滤波得到的xL和xH可以上下排列，因此我们得到了左上角为LL（原始图片的低频近似），左下角为LH（原始图片的水平边），右上角为HL（原始图片的垂直边）和右下角的HH（原始图片的斜边）。
考虑到JPEG压缩算法倾向于保留低频信号、丢弃高频信号，很明显，我们要对LL矩阵动手脚。因为在水平和竖直方向都进行了DWT，因此LL矩阵的长和宽都是原始图片的一半。也正因如此，LL矩阵的4x4像素块就对应了原始图片的8x8像素块。
DCT
离散余弦变换是余弦变换的离散形式（还是废话）。所谓余弦变换，其实就是用余弦来表示原始信号。听起来耳熟不？想想看之前的傅里叶变换，它使用复数表示一个频率的正弦分量（虚部）和余弦分量（实部）。如果简单地说，余弦变换其实就是只保留实部的傅里叶变换。但是从公式上来看，你会发现去掉虚部的傅里叶变换和JPEG压缩使用的DCT还不太一样，主要在于系数。JPEG使用的离散余弦变换得到的结果是正交的。
你要问为什么要正交才行，这我也不知道。但从经验来说，通过某种努力得到的一种特殊性总能带来一些好处。
这里我们使用原始的Orthogonal DCT-II和Orthogonal DCT-III，前者作为DCT，后者作为IDCT。同样地，这是一维变换，但和傅里叶变换类似，DCT也具有分离性，也就是说，我们可以先按行计算DCT，然后将计算结果按列再计算DCT。此外，由于我们只是在4x4的像素块上计算DCT，因此不使用快速算法也能很快地完成运算（使用并行计算）。
当然，使用快速算法固然更好一些。但是只要性能还说得过去，就没有必要优化。
通过对DWT结果的LL矩阵进行DCT运算，我们能够将近似矩阵中的能量分离。虽然LL矩阵本身已经是低频信号了，但是其展现形式仍然是原始图像的近似。通过DCT，我们能够将4x4像素块中的能量以余弦为模板分解。在分解结果中，左上角是低频，右下角是高频。这种位置上的固定能够给后续操作带来很大的便利——我们显然想要在低频中动手脚，因为这是最容易被JPEG压缩算法保留的。
SVD
奇异值分解是一个比较陌生的概念。虽然我很喜欢教我线性代数的老师（北方工业大学的钱盛老师），但很遗憾，因为学完了一直没怎么用到相关的知识，我现在已经完全忘光了。
在开始说奇异值分解之前，我想用特征值来类比一下：方阵A表示一种空间变换，v是向量，$Av=\lambda v$翻译成人话就是向量v经过一种变换A，得到的结果和这个向量经过缩放$\lambda$倍得到的结果相同，在这种情况下，将v视为特征向量，在这个特征向量下$\lambda$就是A的特征值。直观地来说，就是对于一个特征向量，空间变换A和线性缩放$\lambda$的作用效果是一样的。这里$\lambda$描述的是A的一种特性或特征，因此$\lambda$是A的特征值。
说到这里，我记得我肯定是学过求特征值，但现在我真是一步都不记得了。
如果把特征值的定义改写一下：$A = v \lambda v^{T}$，你会发现这个和奇异值分解就很像了：$A = \mu \Sigma \sigma^T$。其中$\mu$和$\sigma$是两组正交单位向量，而$\Sigma$是对角阵，表示奇异值。和特征值的缩放不同，奇异值分解包含了旋转和缩放，其中奇异值就代表了缩放，而余下的两组基则代表了旋转。（可能还有投影，但是考虑到A是一个方阵，所以因该没有投影的事儿）
通过对DCT系数进行奇异值分解，我们能够更进一步地使用奇异值来描述4x4像素块的频率信息。基于能量（频率）分解的角度，我们还可以这样理解奇异值分解：对于$A = U \Sigma V^T$，矩阵U和V的列向量表示一种模式，而$\Sigma$矩阵则描述了这两种模式如何叠加。奇异值矩阵的一些特点使得我们在隐写式更好处理：
- 奇异值按照非递增顺序排列，即从左上角到右下角的方向
- 奇异值是非负的实数，表示矩阵在每个模式或方向上的能量强度。较大的奇异值对应了矩阵中重要的能量模式，而较小的奇异值对应了较不重要的能量模式
基于这两点，再考虑DCT结果近似是一个上三角阵，因此通过奇异值分解后，奇异值矩阵左上角的元素就是值最大的低频信息，也是这个4x4像素块中占主导地位的信息。通过修改这个元素，我们能够将信息隐写到图像的主要信息中，保证了隐写结果能够抵抗JPEG压缩。
编码
虽然我们理清了计算思路，但还有一个问题：如何将bit编码到实数范围的奇异值中？
在上一篇DFT隐写中，我才用了取0和置1的方法——bit0对应频率强度为0；bit1对应频率强度大于一个给定的阈值。这是因为DFT计算结果的差距比较大，且不论正负，在绝对值上就从零点几跨越到了好几万。但这次不一样了——由于DCT和SVD运算只发生在4x4的像素块上，因此元素之间的差距不太大——最大不超过几千。因此我是用余数的方式来编码信息。考虑LSB编码，虽然是直接操作了bit，但换一种思路想，可以认为LSB是取2的余数，如果余数为0，那就是bit0，如果余数为1，那就是bit1。对于奇异值是小数的情况，我们可以选定一个值d，如果是bit 0，那么就让奇异值变成(floor(s / d) + 0.25) * d，这里考虑到SVD、DCT和DWT运算都是小数，计算结果可能会因为精度损失而上下浮动，因此加上0.25让它保证奇异值对d的余数不超过d的一半，同时也不会因为精度损失缩小被反转为bit1 （例如经过取整后应该是20，但经过计算后变成了19.9，原本20对5的余数是0，因此是bit0，但19.9对5的余数就是4.9，变成了bit1）；如果是bit1，则使用(floor(s / d) + 0.75) * d。
在解码时可以直接使用(s % d) / d，如果值大于0.5就可以认定为bit 1，如果值小于0.5就认定为bit 0。当然，如果想让阈值更加灵活，可以使用一维K平均算法找到bit 0和bit 1的中心值，一次来判断是bit 0还是bit 1。
代码实现
代码实现的部分只挑主要的说，完整代码可以移步GitHub仓库。
DWT
我在Java上没找到什么可靠的DWT计算库。有个JWave，但是它的DWT只支持2次幂大小的输入，而我想要支持任意大小的输入。所以我参考了WaveLib，把C代码用Kotlin重写了。
重写的过程肯定不算愉快，但我发现遵循如下技巧能够给自己省不少心：
- 首先不要想C代码是做什么的，原封不动的把代码用Kotlin改写，保证调用起来的效果是一致的。
- 使用调试来确保变量的值是相符的
- 不要试图“融入”语言的特性——人家C怎么写的，你就用别的语言怎么翻译
这里我是用CLion作为IDE调试C代码，因为我对C不是很熟，所以CLion给我提供了很大的帮助，包括根据类型寻找定义，使用调试器打断点和查看变量的值。经过这样一番折腾，我得到了一个丑陋但是和C语言行为一致的代码。
之后我便开始阅读代码，将C代码逐步改写成Kotlin代码，并且保证改写前后代码的行为一致。但最终我发现C语言和Java还是有本质的区别。我注意到C语言中处理二维数组可以像处理一维数组那样。但是在Java中不行——Java中的数组索引只能是Int，这就意味着Java的数组是受到逻辑限制的。因此在Java中该用二维数组就得用二维数组，但是这套C代码是针对一维数组写的，所以在改写的时候花费了不少精力。在C中，通过控制读取数组的步长可以控制DWT是按行（步长为1）还是按列（步长为宽度），在Kotlin中我是用Lambda来实现类似的功能——按行的地方就使用arr[y][it]来访问，按列的地方就是用arr[it][x]来访问。类似地，在C代码中作者复用了数组来存储结果，在理解DWT变换之后，我选择每次计算都产生一个新数组——我都用Java了，我还在乎那点内存和性能吗？
最总我还是成功地把对称填充的DWT2D移植到了Kotlin中，并且使用Kotlin协程加速计算。
由于DWT运算要求二次幂大小的输入，因此需要对原始信号进行填充。这里我选择对称填充。如果原始信号是x1...xn，那么填充后的信号就是...|x3|x2|x1|x1|x2|x3|...|xn|xn|xn-1|...DCT
DCT相对好办，因为尺寸是固定的，所以我只需要按照维基百科上给出的公式计算就行了。其中我注意到一些值其实是可以缓存的，因此我设计了一个查找表。给定一个大小之后（比如4），提前计算出一些三角函数的常量，比起实时运算，很显然直接查表要快非常多。
SVD
之前我使用了Apache Common Math 3的傅里叶计算，我惊喜地发现它也有SVD，因此我就直接用这个库的了。
组装
在有了以上的基础之后，我们就可以把程序组装起来了。程序的流程大概如下：
读入图片，从RGB转换成YUV，如果边长是奇数，则补0变成偶数
- 转换成YUV时，Y通道需要减128使其范围对称，这是DCT的要求
- 分别对YUV三个通道应用Haar DWT，取出LL矩阵，余下三个存起来
- 针对三个LL矩阵，切分成4x4的像素块，不够的地方补0
- 对这些像素块应用DCT和SVD
- 隐写
- 将修改后的S矩阵还原回DCT系数，然后进行IDCT
- 将4x4像素块展开，还原回LL矩阵
- 配合之前存起来的矩阵，使用IDWT还原出隐写后的通道
- 将YUV转换成RGB，保存成图片
在隐写时需要选定参数，但是这个参数会取决于图片本身的特性，甚至不同通道取用的d值都不一样，需要不断尝试。因此我将前四步抽取出来，一张图片计算出SVD之后，隐写时直接给出S矩阵的副本，然后进行后续操作。这样在程序中搜索d值的时候可以免去重复计算一样的东西：
package info.skyblond.steganography.dwtdctsvd
import info.skyblond.steganography.yuvToRGB
import kotlinx.coroutines.async
import kotlinx.coroutines.awaitAll
import kotlinx.coroutines.coroutineScope
import kotlinx.coroutines.runBlocking
import org.apache.commons.math3.linear.RealMatrix
import org.apache.commons.math3.linear.SingularValueDecomposition
import java.awt.image.BufferedImage
class DWTDCTSVD(
private val sourceImage: BufferedImage
) {
private val yDecompositions: List<Decomposition>
private val uDecompositions: List<Decomposition>
private val vDecompositions: List<Decomposition>
private val ySVD: Array<Array<SingularValueDecomposition>>
private val uSVD: Array<Array<SingularValueDecomposition>>
private val vSVD: Array<Array<SingularValueDecomposition>>
private val llWidth: Int
private val llHeight: Int
init {
// RGB -> YUV -> Haar DWT
val (yDWT, uDWT, vDWT) = runBlocking {
decomposeImageYUVHaarDWT(sourceImage)
}
val yLL = yDWT.let { yDecompositions = it.second; it.first }
val uLL = uDWT.let { uDecompositions = it.second; it.first }
val vLL = vDWT.let { vDecompositions = it.second; it.first }
// yuv should share the same size
llWidth = yLL.width
llHeight = yLL.height
// LL -> chunk to 4x4
val yChunks = yLL.chunked4x4()
val uChunks = uLL.chunked4x4()
val vChunks = vLL.chunked4x4()
// Chunk -> DCT -> SVD
ySVD = runBlocking { yChunks.applyDCTSVD() }
uSVD = runBlocking { uChunks.applyDCTSVD() }
vSVD = runBlocking { vChunks.applyDCTSVD() }
}
private suspend fun applyISVDIDCT(
svdChannel: Array<Array<SingularValueDecomposition>>,
modifiedS: Array<Array<RealMatrix>>
): Array<Array<Chunk>> = coroutineScope {
svdChannel.mapIndexed { y, row ->
val rowS = modifiedS[y]
async {
row.mapIndexed { x, svd ->
async {
val recovered = svd.u.multiply(rowS[x]).multiply(svd.vt)
// copy to chunk and apply idct
Chunk(Array(4) { y -> DoubleArray(4) { x -> recovered.getEntry(y, x) } })
.idct2d()
}
}.awaitAll().toTypedArray()
}
}.awaitAll().toTypedArray()
}
private fun copyMatrixS(): Triple<Array<Array<RealMatrix>>, Array<Array<RealMatrix>>, Array<Array<RealMatrix>>> {
val yS = ySVD.let {
Array(it.size) { y -> Array(it[y].size) { x -> it[y][x].s.copy() } }
}
val uS = uSVD.let {
Array(it.size) { y -> Array(it[y].size) { x -> it[y][x].s.copy() } }
}
val vS = vSVD.let {
Array(it.size) { y -> Array(it[y].size) { x -> it[y][x].s.copy() } }
}
return Triple(yS, uS, vS)
}
suspend fun steganography(action: (Array<Array<RealMatrix>>, Array<Array<RealMatrix>>, Array<Array<RealMatrix>>) -> Unit): BufferedImage =
coroutineScope {
// each time we create a copy of SVD matrix, thus we can reuse
val (yS, uS, vS) = copyMatrixS()
action(yS, uS, vS)
// then iSVD -> iDCT -> Chunk
val (yLLChanged, uLLChanged, vLLChanged) = listOf(
async { applyISVDIDCT(ySVD, yS).flatten4x4(llWidth, llHeight) },
async { applyISVDIDCT(uSVD, uS).flatten4x4(llWidth, llHeight) },
async { applyISVDIDCT(vSVD, vS).flatten4x4(llWidth, llHeight) },
).awaitAll()
// LL -> YUV Channel
val yIDWT = idwt2d(yLLChanged, yDecompositions, Wavelet.haar)
val uIDWT = idwt2d(uLLChanged, uDecompositions, Wavelet.haar)
val vIDWT = idwt2d(vLLChanged, vDecompositions, Wavelet.haar)
val recoveredImage = BufferedImage(sourceImage.width, sourceImage.height, sourceImage.type)
for (y in 0 until recoveredImage.height) {
for (x in 0 until recoveredImage.width) {
val cy = yIDWT[y][x] + 128
val cu = uIDWT[y][x]
val cv = vIDWT[y][x]
val color = yuvToRGB(cy, cu, cv)
recoveredImage.setRGB(x, y, color.rgb)
}
}
recoveredImage
}
}由于三个通道的d值并不具有相关性，所以可以单独搜索每个通道的d：
package info.skyblond.steganography.dwtdctsvd
import info.skyblond.steganography.writeJPG
import info.skyblond.steganography.writePNG
import kotlinx.coroutines.Dispatchers
import kotlinx.coroutines.async
import kotlinx.coroutines.awaitAll
import kotlinx.coroutines.runBlocking
import org.apache.commons.math3.linear.RealMatrix
import java.awt.Color
import java.awt.image.BufferedImage
import java.io.File
import javax.imageio.ImageIO
import kotlin.math.abs
import kotlin.math.floor
object Test {
private val sourceImage = ImageIO.read(File("./pic/pexels-lukas-hartmann-1497306.jpg"))
private val outputDir = File("./output/dwt-dct-svd").apply { mkdirs() }
@JvmStatic
fun main(args: Array<String>): Unit = runBlocking(Dispatchers.Default) {
println("Mask: w:${maskImage.width} h:${maskImage.height}")
println("Preparing SVD...")
val obj = DWTDCTSVD(sourceImage)
(20..150 step 2).map { d ->
val encodedImage = searchSingleEncode(obj, d) { y, u, v -> v }
async {
val file = File(outputDir, "test_dY_${dY}_dU_${dU}_dV_$d.jpg")
encodedImage.writeJPG(file, 0.7)
searchSingleDecode(
"jpg_0.7_dY_${dY}_dU_${dU}_dV_$d",
ImageIO.read(file), d
) { y, u, v -> v }
}
}.awaitAll()
}
private suspend fun searchSingleEncode(
obj: DWTDCTSVD, d: Int,
selector: (Array<Array<RealMatrix>>, Array<Array<RealMatrix>>, Array<Array<RealMatrix>>) -> Array<Array<RealMatrix>>
): BufferedImage {
println("Applying steganography... $d")
val iter = mask.iterator()
return obj.steganography { ySArr, uSArr, vSArr ->
for (y in 0 until ySArr.size - 1) {
for (x in 0 until ySArr[0].size - 1) {
val bit = iter.next()
selector(ySArr, uSArr, vSArr)[y][x].setBit(bit, d)
}
}
}
}
private suspend fun searchSingleDecode(
name: String, image: BufferedImage, d: Int,
selector: (SubRegion, SubRegion, SubRegion) -> SubRegion
) {
println("Decoding image... $d")
val (yDWT, uDWT, vDWT) = decomposeImageYUVHaarDWT(image)
val (yLL, _) = yDWT
val (uLL, _) = uDWT
val (vLL, _) = vDWT
val selectedLL = selector(yLL, uLL, vLL)
val selectedDCTSVDChunks = selectedLL.chunked4x4().applyDCTSVD()
val yBitStream = buildList {
for (y in 0 until selectedDCTSVDChunks.size - 1) {
for (x in 0 until selectedDCTSVDChunks[0].size - 1) {
val v = selectedDCTSVDChunks[y][x].s.readBit(d)
add(v)
}
}
}
val data = kMeans1D(yBitStream.toDoubleArray())
val outputImage = BufferedImage(maskImage.width, data.size / maskImage.width, BufferedImage.TYPE_INT_ARGB)
for (y in 0 until outputImage.height) {
for (x in 0 until outputImage.width) {
val color = if (data[y * outputImage.width + x]) Color.WHITE else Color.BLACK
outputImage.setRGB(x, y, color.rgb)
}
}
outputImage.writePNG(File(outputDir, "recovered_$name.png"))
}
/**
* Encode:
* (floor([RealMatrix].getEntry(i,i) / [d]) + f([bit])) * [d],
*
* The result is the SVD entry will be quantized by [d].
* If the bit is true, then the reminder will be bigger than 0.5d_i,
* otherwise it's smaller than 0.5d_i.
* */
private fun RealMatrix.setBit(bit: Boolean, d: Int) {
val b = if (bit) 0.75 else 0.25
val s = this.getEntry(0, 0)
val m = floor(s / d) + b
this.setEntry(0, 0, m * d)
}
private fun RealMatrix.readBit(d: Int): Double {
val s = this.getEntry(0, 0)
return (s % d) / d
}
private fun kMeans1D(inputs: DoubleArray): BooleanArray {
var threshold: Double
val err = 1e-6
// center point of 0 and 1
val center = doubleArrayOf(inputs.min(), inputs.max())
while (true) {
threshold = (center[0] + center[1]) / 2
val bits = BooleanArray(inputs.size) { inputs[it] >= threshold }
// check the distance and correct center
center[0] = inputs.filterIndexed { index, _ -> !bits[index] }.average() // the avg of all bit0's value
center[1] = inputs.filterIndexed { index, _ -> bits[index] }.average() // the avg of all bit1's value
val newThreshold = (center[0] + center[1]) / 2
if (abs(newThreshold - threshold) < err) {
// update the threshold and exit
threshold = newThreshold
break
}
}
return BooleanArray(inputs.size) { inputs[it] >= threshold }
}
private val maskImage = ImageIO.read(File("./pic/fft_mask/mask3.png"))
private val maskData = BooleanArray(maskImage.height * maskImage.width) { index ->
maskImage.getRGB(
index % maskImage.width, index / maskImage.width
) == Color.WHITE.rgb
}
private val mask = object : Iterable<Boolean> {
override fun iterator(): Iterator<Boolean> = object : Iterator<Boolean> {
override fun hasNext(): Boolean = true
private var d = maskData.iterator()
override fun next(): Boolean {
if (!d.hasNext()) d = maskData.iterator()
return d.next()
}
}
}
}这里通过调整searchSingleEncode和searchSingleDecode的最后一个参数就可以选择搜索哪个通道。等搜索到了每个通道表现最佳的d值之后，便可以使用这个函数编解码：
private suspend fun applySteganography(
obj: DWTDCTSVD,
dY: Int, dU: Int, dV: Int,
): BufferedImage {
println("Applying steganography... $dY $dU $dV")
val iter = mask.iterator()
return obj.steganography { ySArr, uSArr, vSArr ->
for (y in 0 until ySArr.size - 1) {
for (x in 0 until ySArr[0].size - 1) {
val bit = iter.next()
uSArr[y][x].setBit(bit, dU)
vSArr[y][x].setBit(bit, dV)
ySArr[y][x].setBit(bit, dY)
}
}
}
}
private suspend fun decode(
name: String, image: BufferedImage,
dY: Int, dU: Int, dV: Int,
) {
println("Decoding image... $dY $dU $dV")
val (yDWT, uDWT, vDWT) = decomposeImageYUVHaarDWT(image)
val (yLL, _) = yDWT
val (uLL, _) = uDWT
val (vLL, _) = vDWT
val yDCTSVDChunks = yLL.chunked4x4().applyDCTSVD()
val uDCTSVDChunks = uLL.chunked4x4().applyDCTSVD()
val vDCTSVDChunks = vLL.chunked4x4().applyDCTSVD()
val yBitStream = buildList {
for (y in 0 until yDCTSVDChunks.size - 1) {
for (x in 0 until yDCTSVDChunks[0].size - 1) {
val vu = uDCTSVDChunks[y][x].s.readBit(dU)
val vv = vDCTSVDChunks[y][x].s.readBit(dV)
val vy = yDCTSVDChunks[y][x].s.readBit(dY)
add(vy + vu + vv)
}
}
}
val data = kMeans1D(yBitStream.toDoubleArray())
val outputImage = BufferedImage(maskImage.width, data.size / maskImage.width, BufferedImage.TYPE_INT_ARGB)
for (y in 0 until outputImage.height) {
for (x in 0 until outputImage.width) {
val color = if (data[y * outputImage.width + x]) Color.WHITE else Color.BLACK
outputImage.setRGB(x, y, color.rgb)
}
}
outputImage.writePNG(File(outputDir, "recovered_$name.png"))
}这里我同时在三个通道进行隐写，然后取平均值。当然你可以使用随机数发生器来决定在哪个通道隐写。不过需要注意的是这并不是一种稳定可靠的隐写方式，因此直接写入字节数据可能会导致数据不可读。
隐写效果
这是一张仅在V通道进行隐写的图片，图片本身经过JPEG压缩，质量大约在70%，这是根据JPEG压缩后的图片解码出来的结果：
如果你自己看的话，你会发现放大之后会有比较明显的格子分界线，同时有些格子比其他格子红一点点：
这是因为在隐写时，选用d本身就是一种量化。我们首先对奇异值进行整除，然后根据bit来决定是加上0.25还是0.75倍的d，这使得d值过大时会产生明显的量化（对于Y通道尤为明显），而d不够大时又没办法保证隐写数据的稳定性。以d=150在Y通道隐写为例：
虽然隐写结果异常清晰，但可以看到隐写后的图片四周出现了明显的条带。
在观察到奇异值的第一个元素带来的改动太大后，我尝试修改了第二个元素，但第二个元素太弱，没办法通过JPEG压缩算法，于是只好作罢。
后记
通过搜索不同的d值，我发现如果你知道你要寻找什么的话，那使用这种隐写方法几乎很难做到修改幅度肉眼不可见，同时还要隐写数据足够稳定。当然，谈及改进的话，我觉得首先可以改进的就是编码方式，很显然不同格子的特性都不一样，因此对于d的选择可以使用自适应算法，根据上下文的相关信息来决定d的取值。不过在设计上有很大一部分受限于无法使用原图。原图的存在会很容易地暴露隐写数据的存在，而不使用原图则要求所有解码用的数据都必须包含在图片里。
我在一开头就说了，这篇文章唯一的依据就是一篇论文的摘要。除此之外全都是我和ChatGPT（主要是我，ChatGPT在冷门领域完全就是个废物，只会复读）填补出来的。经过反复尝试，我发现这种隐写方法并没有达成我的期望——它固然是能够抵抗JPEG压缩，但任何隐写算法只要能够留下足够强的标记，它们总能被JPEG压缩留下。
历时两个月，我想至此关于隐写术就暂且告一段落了。
-全文完-
【代码札记】初探图片隐写术（Steganography）之DWT-DCT-SVD隐写 由 天空 Blond 采用 知识共享 署名 - 非商业性使用 - 相同方式共享 4.0 国际 许可协议进行许可。
本许可协议授权之外的使用权限可以从 https://skyblond.info/about.html 处获得。
论文全文已发送至邮箱。
收到了，感谢。我发现我的方法和论文里说的方法完全不一样（
看不懂，只能膜拜大佬~@(呵呵)