目录

一、背景

二、原理

1、离散Laplace算子介绍

2、Laplace卷积

3、Possion方程解法介绍

 三、验证

 四、Python下的算法实现

a、DCT求解

 1、定义函数calMSE计算误差Mean Square Error

2、导入原图并记下大小

3、拓展原图

4、卷积并求DCT变换

5、除以分母

 6、逆变换、拉回低频并显示。

7、打印误差与耗时

 b、DST求解        

五、误差与结果

六、源码

一、背景

        数字化之后的图像数据占用的空间非常大，我们不妨计算以下，一张分辨率为800*600像素的32位灰度图像，其像素数目为480000，占用空间为480000*32bit=480000*4B≈1.83MB。电影帧率假设24帧每秒，则此分辨率下存储一秒电影需要1.83*24=43.92MB，而目前大多数使用情况确是三通道下的彩图，那么这个数字乘上3，仅仅一秒，可想数据量多么庞大。

        图像数据大小给图像存储及传输制造了很大的困难。图像压缩是将图像存在的冗余信息除去，实现有效压缩。如特征图可以保留各个位置的相对情况，离散傅里叶变换可以把数据集中在低频等等方法。本文仅采用特征图压缩“图像”数据量，利用泊松方程的求解重构图像

二、原理

1、离散Laplace算子介绍

        离散函数的导数退化为差分，一维一阶差分公式和二阶公式分别为：

        分别对Laplace算子x, y两个方向的二阶导数进行差分就得到的离散函数的Laplace算子

而在一个二维函数f（x, y）中，x, y两个方向的二阶差分分别为：

        所以Laplace算子的差分形式为：

         写成卷积核的形式就是：

         当然也可以写成这样：

        注意该核的特点，mask在上下左右四个90°的方向上结果相同，也就是说在90度方向上无方向性。为了让该mask在45°也有该性质，我们引进另一种卷积核：

         这是最常见的两种拉普拉斯卷积核，以下再引进另外一种文中会用到的常见卷积核：

        而所谓卷积核，其操作无非是与空间滤波一样在原图上逐行移动，将核上的数值与其重合的像素相乘后求和，最后将计得的数值赋给核中心所对应的副本位置。而对图像的第一，和最后一行/列，我们无法对其进行统一计算，这时介绍两种解决方案：

1. Laplace图比原图小两行两列或存放0
2. 将原图拓展成M+2，N+2大小的图，等于将原图加一像素的边框，这个框赋值0或赋值最近的行/列

        对应于2中的两种方式，也即人为给定了第一类边界条件（迪利克雷边界）与第二类纽曼边界的条件，这时对应有两种求解。而1解决方案中虽解决了图的大小问题，却未给定任何边界条件，无法通过离散余弦/正弦变换进行数值求解，但可以通过求解系数矩阵的方式求解，本实验利用离散傅里叶变换故而采用第2种方案

2、Laplace卷积

        回归正题，现在给定一个4*4大小的图：

        根据上述介绍的第一种卷积：（原图为OriginalPic以下简记为O）

         我们可以得到以下Laplace图

         这样可以列得4个方程：

        我们要求的是原图各点像素值，但四个方程求解十六个未知数是不可能的，于是我们添加约束条件，比如我们知道十二个边界的值（迪利克雷边界条件），就可以再得12个方程，这样16个方程求解十六个未知数就可以很轻松的利用python求解了。但我们要使用的，是下列介绍的dct与dst变换：

        我们以dct为例，假设已经知道纽曼边界条件，即边界的梯度，我们需要将上述4*4拓展为6*6，如下：

         这样我们求得的特征图就是：

3、Possion方程解法介绍

nullhttps://elonen.iki.fi/code/misc-notes/neumann-cosine/index.html       摘自↑↑↑

        如何通过特征图重建原图？这就要用到泊松方程解法：

        连续域的2D泊松方程形式如下：

         离散域形式为：（μ原图，ρ特征图即上文提到的LaplacePic）

         函数u（x, y）可以表示为：

         于是可以得到：

         我们的目标是求解u。引入两种缩写方便计算：

         列出2D DST与其逆变换2D IDST变换计算公式：

         代入有限差分方程两侧并消除求和符号，得到频域方程：

          又

         所以：

          又

        为了求解μ^，简化方程形式：

         因此：

         代入最初的缩写可以得到：

        根据上述推导，求解泊松方程，总体分为四大步骤：

1. 按照上述方法将原图卷积形成特征图ρ
2. 对ρ求2D DCT得到ρ^
3. 对ρ^每一个像素点除以分母2(cos(Πm/J)+cos(Πn/L)-2)
4. 最后通过对ρ^求2D IDCT变换重建出μ（原图）

 三、验证

         好了,原理已经清楚了，我们来手动验证一下，精度0.1，对上述特征图进行DCT变换得到：   

          除以分母得到：

           2D IDCT变换：

                ​​​​​​​  ​​​​​​​              

        最后别忘了，要减去本图的平均值并加上原图的平均值，这一步的目的是因为该方案只能还原相对量而不能还原低频分量，所以必须要有最后一步拉回水准线的操作。本图平均值为0.65，原图平均值为6.25，拉完并整形化以后为：

          ​​​      ​​​​​​​            

         通过以上求得结果来看，虽然与原图有所出入，但总体相似，这并不能说明该求解方法不好，数据量小，手算误差等都是存在的误差因素，在基数为255的图上，这个差距并不大，下面进行代码实现验证。

 四、Python下的算法实现

a、DCT求解

 1、定义函数calMSE计算误差Mean Square Error

2、导入原图并记下大小

3、拓展原图

4、卷积并求DCT变换

5、除以分母

 6、逆变换、拉回低频并显示。

        这里有一个细节问题，需要把拉高以后的输出图中大于255的赋值为255，小于0的赋值0，否则会出现‘黑极生白，白极生黑’的现象

7、打印误差与耗时

 b、DST求解        

        如此代码实现部分便陈述完毕，不过以上是DCT求解法，对于DST求解，由于相似度极高，只做简单论述。只需将步骤3改为：

         即将拓展的边界赋值0，利用迪利克雷边界条件求解。再将步骤5改为：

        即对该条件下的离散正弦求解。最后将DCT中用到的两次正逆变换改为DST正逆变换即可。但是我从cv2中并没有找到dst()函数，于是我们使用scipy库里的fft.dstn()与fft.idstn()函数替换。

        上述两种边界条件下的求解方法虽然都可以对特征图进行还原，但还原精确度却大不相同，实验验证，在不同的图像下两种重建方法处理速度基本相同。

        至此，整个实验步骤结束

五、误差与结果

     DCT下：

        原图：

        ​​​​​​​        ​​​​​​​        

特征图：

        ​​​​​​​        ​​​​​​​        

第一种卷积核：（MSE三个值分别对应三个通道的均方差）

        ​​​​​​​        ​​​​​​​        

 耗时与误差：

​​​​​​​​​​​​​​

均方差控制在0.4以内，图像还原很成功。

第二种：（颜色很暗，还原度不高）

        ​​​​​​​        ​​​​​​​        

 耗时与误差：

 第三种：

        ​​​​​​​        ​​​​​​​        

 误差与耗时：

 六、源码

import numpy as np
import cv2
import time
from scipy import fft

time_start = time.time()  # 开始计时


# 定义函数calMSE，计恢复图std与目标图像aim之间的均方误差（MSE）评价指标
def calMSE(std, aim):
    # 循环计算每个通道均方差
    num = np.zeros((std.shape[2]), dtype=np.float32)
    mse = np.zeros((std.shape[2]), dtype=np.float32)
    for c in range(std.shape[2]):
        for i in range(std.shape[0]):
            for j in range(std.shape[1]):
                num[c] += np.power(np.float32(std[i, j, c]) - np.float32(aim[i, j, c]), 2)  # 这里需要强制转换
        mse[c] = np.sqrt(num[c] / (std.shape[0] * std.shape[1]))  # 否则误差偏小
    return mse


def main():
    img_cv = cv2.imread("D:/sy15.jpg")

    OriginalPic0 = np.array(img_cv, dtype=np.uint8)
    M, N, ch = OriginalPic0.shape  # M,N--长、宽（单通道）

    OriginalPic = np.zeros((OriginalPic0.shape[0] + 2, OriginalPic0.shape[1] + 2, 3), dtype=np.uint8)
    img_dct = np.zeros((M, N, ch), dtype=np.float32)
    img_idct = np.zeros((M, N, ch), dtype=np.float32)

    for c in range(ch):
        for i in range(1, M + 1):
            for j in range(1, N + 1):
                OriginalPic[i, j, c] = OriginalPic0[i - 1, j - 1, c]
        OriginalPic[:, 0, c] = OriginalPic[:, 1, c]
        OriginalPic[:, N + 1, c] = OriginalPic[:, N, c]
        OriginalPic[0, :, c] = OriginalPic[1, :, c]
        OriginalPic[M + 1, :, c] = OriginalPic[M, :, c]

    LaplacePic = np.zeros((M, N, ch), dtype=np.float32)  # 存放特征图
    LaplacePic_show = np.zeros((M, N, ch), dtype=np.uint8)  # 待显示的特征图

    img_recor = np.zeros((M, N, ch), dtype=np.float32)

    OriginalMeanValue = np.zeros(3, dtype=np.float32)
    for c in range(ch):
        OriginalMeanValue[c] = np.sum(OriginalPic[:, :, c]) / (M * N)
    print("OriginalMeanValue：", OriginalMeanValue)

    kernel1 = [[0, 1, 0], [1, -4, 1], [0, 1, 0]]  # 拉普拉斯几种卷积
    kernel2 = [[1 / 4, 1 / 4, 1 / 4], [1 / 4, -2, 1 / 4], [1 / 4, 1 / 4, 1 / 4]]  # 这里因为还原后色差不明显
    kernel3 = [[1 / 4, 1 / 2, 1 / 4], [1 / 2, -3, 1 / 2], [1 / 4, 1 / 2, 1 / 4]]  # 手动按比例缩放

    for c in range(ch):
        for i in range(1, M + 1):
            for j in range(1, N + 1):
                LaplacePic[i - 1][j - 1][c] = (np.sum(np.multiply(kernel3, OriginalPic[i - 1:i + 2, j - 1:j + 2, c])))
                LaplacePic_show[i - 1][j - 1][c] = LaplacePic[i - 1][j - 1][c] + 128  # 特征明显化
        img_dct[:, :, c] = cv2.dct(LaplacePic[:, :, c])
    cv2.imshow("Original", OriginalPic)  # 显示原图
    cv2.imshow("Laplace", LaplacePic_show)  # 显示特征图

    print("Original", OriginalPic0.shape)
    print("Laplace", LaplacePic.shape)

    for c in range(ch):
        for i in range(0, M):
            for j in range(0, N):
                img_dct[i, j, c] /= (2 * (np.cos(np.pi * i / (M + 2)) + np.cos(np.pi * j / (N + 2)) - 2))  # 除以2cosΠ...
        img_dct[0, 0, c] = 0
        img_idct[:, :, c] = cv2.idct(img_dct[:, :, c])


    IdctMeanValue = np.zeros(3, dtype=np.float32)
    for c in range(ch):
        IdctMeanValue[c] = np.sum(img_idct[:, :, c]) / (M * N)
    print("IdctMeanValue：", IdctMeanValue)
    for c in range(ch):
        for i in range(0, M):
            for j in range(0, N):
                img_recor[i, j, c] = img_idct[i, j, c] + OriginalMeanValue[c] - IdctMeanValue[c]
                if (img_recor[i, j, c] > 255):
                    img_recor[i, j, c] = 255
                if (img_recor[i, j, c] < 0):
                    img_recor[i, j, c] = 0

    img_recor = np.uint8(img_recor)  # 在这之前要将大于255的定义255，小于零的定为0
    print("recor", img_recor.shape)  # 否则-1 -> 254 原黑变白
    cv2.imshow("img_recor", img_recor)

    MeanSquareError = calMSE(img_recor, OriginalPic0)
    print("MeanSquareError", MeanSquareError)

    time_end = time.time()  # 结束计时
    print("Time cost = %fs" % (time_end - time_start))
    cv2.waitKey(0)


if __name__ == '__main__':
    main()