离散域下的泊松方程求解(python实现)

目录

一、背景

二、原理

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算子介绍

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

                            离散域下的泊松方程求解(python实现)

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

而在一个二维函数f(x, y)中,x, y两个方向的二阶差分分别为:

                          离散域下的泊松方程求解(python实现)

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

     离散域下的泊松方程求解(python实现)

         写成卷积核的形式就是:

                              离散域下的泊松方程求解(python实现)

         当然也可以写成这样:

                                离散域下的泊松方程求解(python实现)

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

                                离散域下的泊松方程求解(python实现)

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

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

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

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

2、Laplace卷积

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

                          离散域下的泊松方程求解(python实现)

        根据上述介绍的第一种卷积:(原图为OriginalPic以下简记为O)

              离散域下的泊松方程求解(python实现)

         我们可以得到以下Laplace图

                                            离散域下的泊松方程求解(python实现)

         这样可以列得4个方程:

                    离散域下的泊松方程求解(python实现)

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

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

          离散域下的泊松方程求解(python实现)

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

                  离散域下的泊松方程求解(python实现)

 

3、Possion方程解法介绍

null离散域下的泊松方程求解(python实现)https://elonen.iki.fi/code/misc-notes/neumann-cosine/index.html       摘自↑↑↑

        

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

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

                                                  离散域下的泊松方程求解(python实现)

         离散域形式为:(μ原图,ρ特征图即上文提到的LaplacePic)

                      离散域下的泊松方程求解(python实现)

         函数u(x, y)可以表示为:

                                  离散域下的泊松方程求解(python实现)

         于是可以得到:

              离散域下的泊松方程求解(python实现)

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

                                                         离散域下的泊松方程求解(python实现)

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

                                                     离散域下的泊松方程求解(python实现)

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

离散域下的泊松方程求解(python实现)

          又

                                 离散域下的泊松方程求解(python实现)

         所以:

                                   离散域下的泊松方程求解(python实现)

          又

                                              离散域下的泊松方程求解(python实现)

        为了求解μ^,简化方程形式:

            离散域下的泊松方程求解(python实现)

             离散域下的泊松方程求解(python实现)

         因此:

      离散域下的泊松方程求解(python实现)

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

                            离散域下的泊松方程求解(python实现)

 

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

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

 三、验证

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

                                     离散域下的泊松方程求解(python实现)

          除以分母得到:

                                离散域下的泊松方程求解(python实现)

           2D IDCT变换:

                ​​​​​​​  ​​​​​​​              离散域下的泊松方程求解(python实现)

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

          ​​​      ​​​​​​​            离散域下的泊松方程求解(python实现)

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

 

 四、Python下的算法实现

a、DCT求解

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

离散域下的泊松方程求解(python实现)

 

2、导入原图并记下大小

离散域下的泊松方程求解(python实现)

 

3、拓展原图

离散域下的泊松方程求解(python实现)

4、卷积并求DCT变换

离散域下的泊松方程求解(python实现)

5、除以分母

离散域下的泊松方程求解(python实现)

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

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

离散域下的泊松方程求解(python实现)

7、打印误差与耗时

离散域下的泊松方程求解(python实现)

 b、DST求解        

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

离散域下的泊松方程求解(python实现)

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

离散域下的泊松方程求解(python实现)

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

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

        至此,整个实验步骤结束

五、误差与结果

     DCT下:

        原图:

        ​​​​​​​        ​​​​​​​        离散域下的泊松方程求解(python实现)

特征图:

        ​​​​​​​        ​​​​​​​        离散域下的泊松方程求解(python实现)

第一种卷积核:(MSE三个值分别对应三个通道的均方差)

        ​​​​​​​        ​​​​​​​        离散域下的泊松方程求解(python实现)

 耗时与误差:

​​​​​​​​​​​​​​离散域下的泊松方程求解(python实现)

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

第二种:(颜色很暗,还原度不高)

        ​​​​​​​        ​​​​​​​        离散域下的泊松方程求解(python实现)

 耗时与误差:离散域下的泊松方程求解(python实现)

 第三种:

        ​​​​​​​        ​​​​​​​        离散域下的泊松方程求解(python实现)

 误差与耗时:

离散域下的泊松方程求解(python实现)

 六、源码

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()

上一篇:numpy实现数学中的各种积


下一篇:numpy 数据类型和空值