盲去模糊blind deconv

元气少女缘结神

已于 2023-03-01 18:01:12 修改

阅读量1.2k

点赞数

分类专栏： Image Deblur 文章标签： opencv

于 2015-04-25 17:22:33 首次发布

本文链接：https://blog.csdn.net/wd1603926823/article/details/45271889

版权

Image Deblur 专栏收录该内容

3 篇文章 0 订阅

订阅专栏

基本概念

看了一些盲去模糊算法，一般主要分为两个步骤：一是从模糊图中计算出核，然后再用这个核去模糊图处理得到清晰图结果。
参考：https://blog.csdn.net/qq_42604176/article/details/105160993
1，对一幅图img进行傅里叶变换后可得到其频域图像，而此过程一般是dft，速度慢；为了提速，一般对原图边界进行扩充预处理，使得size更适合快速傅里叶变换（有现成的第三方库 fftw3，linux下可直接安装）；
2，无论是哪种傅里叶变换，都只能处理浮点型数据；
3，使用图像kernel对图像img进行卷积就等同于：对kernel进行傅里叶变换kernel_dft或kernel_fft，对img进行傅里叶变换img_dft或img_fft，然后将两者傅里叶变换结果点积（常使用此函数cv::mulSpectrums）。

dft-idft小实例

以下的小实例展示了对图像进行dft、显示频谱图、再idft、再显示离散傅里叶逆变换的结果图的过程：

//参考：https://blog.csdn.net/zhaitianbao/article/details/117984283
#include<iostream>
#include<opencv2/opencv.hpp>
#include<ctime>
using namespace std;
using namespace cv;

void fftshift(cv::Mat &plane0, cv::Mat &plane1);

int main()
{
	Mat test = imread("G:\\projects\\XRT\\images\\sourceimgs\\deblur\\0.png", 0);
	test.convertTo(test, CV_32FC1);

	//创建通道，存储dft后的实部与虚部（CV_32F，必须为单通道数）
	cv::Mat plane[] = { test.clone(), cv::Mat::zeros(test.size() , CV_32FC1) };
	cv::Mat complexIm;
	cv::merge(plane, 2, complexIm); // 合并通道 （把两个矩阵合并为一个2通道的Mat类容器）
	cv::dft(complexIm, complexIm, 0); // 进行傅立叶变换，结果保存在自身

    // 分离通道（数组分离）
	cv::split(complexIm, plane);
	// 以下的操作是频域迁移 将图像的高频分量平移到图像的中心
	fftshift(plane[0], plane[1]);
	// 计算幅值
	cv::Mat mag,mag_log,mag_nor;
	cv::magnitude(plane[0], plane[1], mag);
	// 幅值对数化：log（1+m），便于观察频谱信息
	mag += Scalar::all(1);
	cv::log(mag, mag_log);
	cv::normalize(mag_log, mag_nor, 1,0, NORM_MINMAX);
	imwrite("fft.png",mag_nor*255);

	// 再次搬移回来进行逆变换
	fftshift(plane[0], plane[1]);
	cv::Mat BLUR;
	cv::merge(plane, 2, BLUR); // 实部与虚部合并
	cv::idft(BLUR, BLUR);       // idft结果也为复数
	BLUR = BLUR / BLUR.rows / BLUR.cols;
	cv::split(BLUR, plane);//分离通道，主要获取通道

	imwrite("original.png", test);
	imwrite("result.png", plane[0]);

	return 0;
}

// fft变换后进行频谱搬移 plane0:实部   plane1：虚部
void fftshift(cv::Mat &plane0, cv::Mat &plane1)
{
	// 以下的操作是移动图像  (零频移到中心)
	int cx = plane0.cols / 2;
	int cy = plane0.rows / 2;
	cv::Mat part1_r(plane0, cv::Rect(0, 0, cx, cy));  // 元素坐标表示为(cx, cy)
	cv::Mat part2_r(plane0, cv::Rect(cx, 0, cx, cy));
	cv::Mat part3_r(plane0, cv::Rect(0, cy, cx, cy));
	cv::Mat part4_r(plane0, cv::Rect(cx, cy, cx, cy));

	cv::Mat temp;
	part1_r.copyTo(temp);  //左上与右下交换位置(实部)
	part4_r.copyTo(part1_r);
	temp.copyTo(part4_r);

	part2_r.copyTo(temp);  //右上与左下交换位置(实部)
	part3_r.copyTo(part2_r);
	temp.copyTo(part3_r);

	cv::Mat part1_i(plane1, cv::Rect(0, 0, cx, cy));  //元素坐标(cx,cy)
	cv::Mat part2_i(plane1, cv::Rect(cx, 0, cx, cy));
	cv::Mat part3_i(plane1, cv::Rect(0, cy, cx, cy));
	cv::Mat part4_i(plane1, cv::Rect(cx, cy, cx, cy));

	part1_i.copyTo(temp);  //左上与右下交换位置(虚部)
	part4_i.copyTo(part1_i);
	temp.copyTo(part4_i);

	part2_i.copyTo(temp);  //右上与左下交换位置(虚部)
	part3_i.copyTo(part2_i);
	temp.copyTo(part3_i);
}

以上是opencv的dft小实例的简短过程（dft的结果做频谱搬移的目的很多都是为了去噪，因为低频代表了图像的很多有用信息，而其原本分布在图像的四个角附近，高频代表噪声位于图像中间，为了很好的显示图像，我们通常将其换个位置，另低频信息在图像中间。同时另高频信息为0，这样就可以简单去噪）。如果不用opencv则可以用第三方库fftwf，以下是fftwf这个库的介绍：
http://www.4k8k.xyz/article/u011389706/106306166 大家可以看看。

简单的去模糊

先将图像转成float类型，然后对每列做如下操作即得到去模糊后的图：

#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/core/core.hpp>

#include <stdint.h>
#include <stdbool.h>

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

using namespace cv;
using namespace std;

void turnToFloat(Mat input_v,Mat &standardimg);
void deconv_RL(Mat &input,float *psfkernel,Size psfsize,int iterations,float gradientdownrate);

float * conv_same(float *array1,int len1,float *array2,int len2);

int main()
{
	char srcimgfolder[400]={0};

	const int psfheight=5;

	Size psf_size(1,psfheight);
	float psf_kernel[psfheight]={0.11419,0.23412,0.30337,0.23412,0.11419};

	int iterations=5;
	float gradientdownrate=2.5f;

	for(int index=1;index<=1;index++)
	{
		sprintf(srcimgfolder,"G:\\projects\\XRT\\images\\sourceimgs\\deblur\\xrt_images\\%d_low.tiff",index);
		Mat blurredimg=imread(srcimgfolder,-1);//0);

		if(!(blurredimg.data))
			continue;

		Mat blurfloat;
		turnToFloat(blurredimg,blurfloat);

		deconv_RL(blurfloat,psf_kernel,psf_size,iterations,gradientdownrate);

		Mat outimg=Mat::zeros(blurredimg.size(),CV_16UC1);
		for(int r=0;r!=outimg.rows;r++)
		{
			for(int c=0;c!=outimg.cols;c++)
			{
				int tmp=65536*(blurfloat.at<float>(r,c));

				outimg.at<ushort>(r,c)=(tmp>52428)?52428:tmp;
			}
		}

		imwrite("deblurred.tiff",outimg);
	}

	return 0;
}


void deconv_RL(Mat &input,float *psfkernel,Size psfsize,int iterations,float gradientdownrate)
{
	int inrows=input.rows;
	int incols=input.cols;

	float *coldata=(float*)malloc(inrows*sizeof(float));

	for(int c=0;c!=incols;c++)
	{
		//对每列数据进行处理
		double maxvalue=0;
		for(int id=0;id!=inrows;id++)
		{
			coldata[id] = input.ptr<float>(id)[c];
			if(coldata[id]>maxvalue)
			{
				maxvalue=coldata[id];
			}
		}

		for(int iter=0;iter!=iterations;iter++)
		{
			float *conv1st=conv_same(coldata,inrows,psfkernel,psfsize.height);

			for(int conv1stid=0;conv1stid!=inrows;conv1stid++)
			{
				conv1st[conv1stid]=(input.ptr<float>(conv1stid)[c]) / conv1st[conv1stid] ;
			}

			float *conv2nd=conv_same(conv1st,inrows,psfkernel,psfsize.height);

			for(int conv1stid=0;conv1stid!=inrows;conv1stid++)
			{
				float tmp=(coldata[conv1stid])*pow(conv2nd[conv1stid],gradientdownrate);
				if(tmp<0)
				{
					tmp=0;
				}
				if(tmp>maxvalue)
				{
					tmp=maxvalue;
				}
				coldata[conv1stid]=tmp;
			}
		}
		//show
		for(int id=0;id!=inrows;id++)
		{
			input.ptr<float>(id)[c]=coldata[id];
		}
	}

}


//array1,array2 are one dimension. len1>=len2
float * conv_same(float *array1,int len1,float *array2,int len2)
{
	int result_len=len1+len2-1;

	float *tmpout=(float*)malloc(result_len*sizeof(float));
	memset(tmpout,0.f,result_len*sizeof(float));

	for(int id1=0;id1!=len1;id1++)
	{
		for(int id2=0;id2!=len2;id2++)
		{
			tmpout[id1+id2] += (array1[id1]*array2[id2]);
		}
	}

	int lenmax=max(len1,len2);
	int lenmin=min(len1,len2);
	float *out=(float*)malloc(lenmax*sizeof(float));
	for(int id=0;id!=lenmax;id++)
	{
		out[id] = tmpout[id+lenmin/2];
	}

	return out;
}


void turnToFloat(Mat input_v,Mat &standardimg)
{
	int rowsh=input_v.rows;
	int colsw=input_v.cols;

	double minval=0,maxval=0;
	cv::minMaxLoc(input_v,&minval,&maxval);

	Mat standardimgv(rowsh,colsw,CV_32FC1);
	for(int r=0;r!=rowsh;r++)
	{
		for(int c=0;c!=colsw;c++)
		{
			standardimgv.at<float>(r,c)=float(input_v.at<ushort>(r,c))/65535;
		}
	}
	standardimgv.copyTo(standardimg);
}

参考：https://blog.csdn.net/qq_45732223/article/details/110754731
conv--same
但这个算法太简单了，固定的核，注定了适用性比较差，对有的图去模糊效果较好，但对有的图却效果较差，如下去模糊的效果图就非常差：
在这里插入图片描述

Split Bregman

由此可知我们需要更精确的模糊核，直接为每一张模糊图计算其特定的模糊核，如下单尺度计算模糊核，然后利用优化的分裂布雷格曼算法去模糊得到最终结果图如下，最上面是模糊图，中间到下面是去模糊程度不同的结果，可以看到最下面的最清晰，而且噪声不大：

下面又用多尺度计算模糊核然后去模糊的结果如下，左边是原图右边是去模糊结果，可以看到当使用多尺度去模糊时，虽然图像清晰了，但引入了许多噪声：

所以对于这种图像来说，可能多尺度去模糊并不合适，单尺度更合适。代码如

https://download.csdn.net/download/wd1603926823/87517539 中所示。

Convolutional Deblurring for Natural Imaging

我看大多数的盲去模糊算法都非常耗时，因为FFT本来就是耗时的操作，哪怕使用fftwf加速库，仍旧很慢，况且我们常常需要多次迭代，这就造成了最终的算法整体性能低。所以最优的方案就是使得傅里叶变换次数尽量少，比如像《Convolution Deblurring for Natural Imaging》中所述，它的步骤是：

一、将输入的模糊图归一化为0~1之间的float类型图，并按0扩充为奇数行、奇数列的待处理图；

二、再对待处理图进行直接下采样，或对其低频图进行下采样得到下采样图片；

三、计算待处理图的频谱、下采样图的频谱，根据下图将频谱从二维变成1维：

比如从0到CV_PI均分成多少份（作为x），每份上都有对应的频谱值（作为y），然后进行高斯拟合得到最终的广义高斯模型GGaussian的各个系数。最后根据系数可算出psf，即估计的模糊核。

四、根据上述得到的模糊核计算去模糊核deblurring_kernel；

五、对去模糊核deblurring_kernel使用GGaussian进行卷积得到去模糊核的逆inverse_deblurring_kernel，即下图的D；

六、对输入的模糊图fB使用inverse_deblurring_kernel进行conv就相当于对其x方向去模糊即下图的fB*Dx；那么将inverse_deblurring_kernel转置后再做一遍相当于得到了fB*Dxy；直接对输入的模糊图使用转置的去模糊核卷积就相当于得到fB*Dy；

七、上图中还差系数gamma，按下图计算这个系数：相当于两个熵的比值，可以用来调节去模糊的程度。

至此得到了去模糊的结果图fR。这个算法的整体性能比之前的都高。

图像增强

使用图像增强/锐化来代替去模糊，达到肉眼上变清晰的效果，如下所示：

可以看到右边比左边清晰了很多，也没有deblur常见的振铃现象。经常使用的就是Laplace锐化算子对模糊图增强，达到去模糊效果。这样最大的优点就是速度快。

Mat sharpen_op = (Mat_<float>(3, 3) << 0, -1, 0,
						-1, 5, -1,
						0, -1, 0);
		filter2D(srcimg, result1, CV_32F, sharpen_op);

元气少女缘结神

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
打赏
1
评论
盲去模糊blind deconv

3，使用图像kernel对图像img进行卷积就等同于：对kernel进行傅里叶变换kernel_dft或kernel_fft，对img进行傅里叶变换img_dft或img_fft，然后将两者傅里叶变换结果点积（常使用此函数cv::mulSpectrums）。参考：https://blog.csdn.net/qq_42604176/article/details/105160993。参考：https://blog.csdn.net/qq_45732223/article/details/110754731。
复制链接

扫一扫