数字化之后的图像数据占用的空间非常大，我们不妨计算以下，一张分辨率为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°也有该性质，我们引进另一种卷积核：

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

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

Laplace图比原图小两行两列或存放0
将原图拓展成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变换计算公式：

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

又

所以：

又

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

因此：

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

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

按照上述方法将原图卷积形成特征图ρ
对ρ求2D DCT得到ρ^
对ρ^每一个像素点除以分母2(cos(Πm/J)+cos(Πn/L)-2)
最后通过对ρ^求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()

标签：原图,泊松,OriginalPic,img,求解,python,离散,range,np
来源： https://blog.csdn.net/qq_57731481/article/details/122370163

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

离散域下的泊松方程求解（python实现）

一、背景