离散域下的泊松方程求解(python实现) | 您所在的位置:网站首页 › python求奇偶数 › 离散域下的泊松方程求解(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算子介绍离散函数的导数退化为差分,一维一阶差分公式和二阶公式分别为: 分别对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方程解法介绍 null
如何通过特征图重建原图?这就要用到泊松方程解法: 连续域的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、拓展原图 这里有一个细节问题,需要把拉高以后的输出图中大于255的赋值为255,小于0的赋值0,否则会出现‘黑极生白,白极生黑’的现象 如此代码实现部分便陈述完毕,不过以上是DCT求解法,对于DST求解,由于相似度极高,只做简单论述。只需将步骤3改为: 即将拓展的边界赋值0,利用迪利克雷边界条件求解。再将步骤5改为: 即对该条件下的离散正弦求解。最后将DCT中用到的两次正逆变换改为DST正逆变换即可。但是我从cv2中并没有找到dst()函数,于是我们使用scipy库里的fft.dstn()与fft.idstn()函数替换。 上述两种边界条件下的求解方法虽然都可以对特征图进行还原,但还原精确度却大不相同,实验验证,在不同的图像下两种重建方法处理速度基本相同。 至此,整个实验步骤结束 五、误差与结果DCT下: 原图: 特征图: 第一种卷积核:(MSE三个值分别对应三个通道的均方差) 耗时与误差: 均方差控制在0.4以内,图像还原很成功。 第二种:(颜色很暗,还原度不高) 耗时与误差: 第三种: 误差与耗时: |
CopyRight 2018-2019 实验室设备网 版权所有 |