2013-04-10 14:51:37 mlkiller 阅读数 4422

本章主要讲图像处理中的模糊处理部分

英文叫做blur, 也叫做smootiing,  中文中叫做模糊或者平滑。

用过photoshop的人都应该知道,滤镜里面就有模糊这个选项,我们现在看看它是怎么实现的。

一含义

   模糊(平滑)是一种常用的图片处理方式,它的作用可以用来降低噪声,还有其他用途

   看一下opencv 里面的公式

               g(i,j) = \sum_{k,l} f(i+k, j+l) h(k,l)

     g(i,j)是目标坐标的像素值, f(i+k,j+l)是k,l这些地方的像素值, h(k,l)是 kernel,  我不知道怎么去准确翻译它的意义,它是过滤器的系数。 

    简单的按照我的思路去理解,就是一个权值,模糊的含义是将所有的像素按照一定的权值进行运算,得到一个比较均衡的结果。

二 类型

类型有很多种:
均值模糊(box blur) 高斯模糊(gaussian blur)  中值模糊(media blur) 二值模糊(bilateral blur)
本文只讲均值模糊和高斯模糊

三 算法

1 均值模糊
   均值模糊很简单就是周边所有的影响都是1,求平均值即可
K = \dfrac{1}{K_{width} \cdot K_{height}} \begin{bmatrix}    1 & 1 & 1 & ... & 1 \\    1 & 1 & 1 & ... & 1 \\    . & . & . & ... & 1 \\    . & . & . & ... & 1 \\    1 & 1 & 1 & ... & 1   \end{bmatrix}
2 高斯模糊
关于高斯模糊的算法,推荐这个文章
根据这个公式计算出系数即可。
上篇文章写得很详细,我就不班门弄斧了。

四均值模糊的代码和效果

     先放上均值模糊的代码
void boxblur(Mat input ,Mat &out, int x, int y)
{
	// accept only char type matrices
	CV_Assert(input.depth() != sizeof(uchar));

	out.create(input.size(),input.type());

	int nChannels = input.channels();
	int nRows = input.rows;
	int nCols = input.cols;

	int size = x * y;
	float kernel = 1.0/size;

	int i,j;
	uchar* p;
	uchar* q;
	uchar R,G,B;

	for( i = x; i < nRows - x; ++i)
	{
		q = out.ptr<uchar>(i);
		for ( j = y; j < nCols - y; ++j)
		{
			float sumR = 0;
			float sumG = 0;
			float sumB = 0;
			for (int k =0; k<x;k++)
			{
				p = input.ptr<uchar>(i-x+k);
				for(int l = 0; l < y;l++)
				{
					sumB += input.at<uchar>(i - x + k,(j + l - y)*nChannels) * kernel;//p[(l + j -y)*nChannels ] * kernel;
					sumG += input.at<uchar>(i - x + k,(j + l - y)*nChannels + 1) * kernel;//p[(l + j -y)*nChannels + 1] * kernel;
					sumR += input.at<uchar>(i - x + k,(j + l - y)*nChannels + 2) * kernel;//p[(l + j -y)*nChannels + 2] * kernel;
				}
			}
			q[j*nChannels] = sumB;
			q[j*nChannels+1] = sumG;
			q[j*nChannels+2] = sumR;
		}
	}


}

红色部分是我想直接用at,而不用指针,但是效率低的厉害。


下图是用指针的相差了20倍。。。可见指针虽然万恶,但是确实是个好东西。



由于size(4,4)图太小看不清, 实际用的是8
原始 opencv 本文


五高斯模糊的代码和效果

代码如下:

void gaussblur(Mat input ,Mat &out, int x, int y)
{
	float sigma = 1.5;
	Mat kernel;
	float pi = 3.1415926;

	kernel.create(x ,y ,CV_32F);

	float mx = x/2.0;
	float my = y/2.0;

       //这里有问题,后面做修正。
	for (int i =0; i< x;i++)
	{
		for (int j =0; j<y;j++)
		{
			kernel.at<float>(i,j) = exp(-1 * ((i - mx) * (i - mx) +(j - my) * (j-my) )/( 2 * sigma * sigma))/(2 * pi * sigma *sigma) ;
		}
	}


    int nChannels = input.channels();
	int nRows = input.rows;
	int nCols = input.cols;

	out.create(input.size(),input.type());
    uchar* p;
	uchar* q;
	float* s;

	for(int  i = x; i < nRows - x; ++i)
	{
		q = out.ptr<uchar>(i);
		for (int j = y; j < nCols - y; ++j)
		{
			float sumR = 0;
			float sumG = 0;
			float sumB = 0;
			for (int k =0; k<x;k++)
			{
				p = input.ptr<uchar>(i-x+k);
				s = kernel.ptr<float>(k); 
				for(int l = 0; l < y;l++)
				{
					sumB += p[(l + j -y)*nChannels ] * s[l];//input.at<uchar>(i - x + k,(j + l - y)*nChannels) * kernel;//
					sumG += p[(l + j -y)*nChannels + 1] *s[l];//input.at<uchar>(i - x + k,(j + l - y)*nChannels + 1) * kernel;//
					sumR += p[(l + j -y)*nChannels + 2] * s[l];//input.at<uchar>(i - x + k,(j + l - y)*nChannels + 2) * kernel;
				}
			}
			q[j*nChannels] = sumB;
			q[j*nChannels+1] = sumG;
			q[j*nChannels+2] = sumR;
		}
	}

	
}

效率如下:

效果图如下:
本文没有考虑边界的情况,所以都是灰色的,可以考虑一下如何处理边界。
原始 opencv 本文

上面代码有两处问题:
第一是在size比较小的时候,这些点的概率之和不等于1,会导致图片出问题。修正如下:

	float sum = 0;
	for (int i =0; i< x;i++)
	{
		for (int j =0; j<y;j++)
		{
			sum+= kernel.at<float>(i,j) = exp(-1 * ((i - mx) * (i - mx) +(j - my) * (j-my) )/( 2 * sigma * sigma))/(2 * pi * sigma *sigma) ;
		}
	}
	for (int i =0; i< x;i++)
	{
		for (int j =0; j<y;j++)
		{
			kernel.at<float>(i,j) = kernel.at<float>(i,j)/ sum ;
		}
	}


第二个问题是本文中sigma 是个固定值,实际上它是个可变值,具体怎么计算,我没有搞清楚,可以查看opencv的源代码,下面文章有参考价值

更新一下参考opencv里面的可以这样计算
sigma = 0.3*((ksize-1)*0.5 - 1) + 0.8 .
修改程序之后发现和原始的高斯函数基本一致,希望广大朋友们多多评论,本人水平有限,很多地方有纰漏,希望能够共同提高。
2019-08-05 17:28:08 z634863434 阅读数 46

前言

本文使用的环境为:Qt5.11 + OpenCV3.4.6
环境安装参考文档:https://blog.csdn.net/z634863434/article/details/89950961

概念

图像模糊从字面上理解,就是将一张清晰的图像变的模糊不清。在图像处理中,模糊可以理解为对每个像素进行滤波或者平滑处理,使得图像内部和边缘都变得平滑,边界不清晰。图像模拟主要可以用来突显出图像中的明显的特征点,通过模糊我们可以对图像进行特征点的提取或者做运动模糊的功能。而在做图像模糊处理时候,其本质为给图像进行降噪处理,在数学上使用的卷积方式实现,常用的有以下4种方法:均值滤波、高斯滤波、中值滤波以及双边滤波。

均值滤波(归一化盒子滤波)

均值滤波取掩模矩阵大小内所有像素点的平均值赋予给中心像素点,该方法对于边缘值无法进行处理。具体解释如下:
假设有6X6的图像像素点矩阵,在6X6矩阵上有一个3X3的小矩阵,依次从左右向右,再从上往下移动。将3X3矩阵中所有点求和取平均值后,赋值给3X3中的红点,向左移动直到将6X6图像中的像素点全部处理完成。
在这里插入图片描述
数学公式如下:
在这里插入图片描述
如上图1中心点数经过均值模糊后为(2+6+4+1+7+3+1+2+2)/ 9 = 3

API函数:

void cv::blur ( InputArray  src,
				OutputArray  dst,
				Size  	ksize,
				Point  	anchor = Point(-1,-1),
				int  	borderType = BORDER_DEFAULT 
	) 	

作用:均值滤波,图像均值模糊处理

输入参数 参数定义
src 输入图像
dst 输出图像
ksize 模糊核大小
anchor 锚点,(-1,-1)为中心点
borderType 边框模式用于推断图像外部的像素

高斯滤波

高斯滤波对图像处理的过程与均值滤波类似,但高斯滤波符合正态分布,在计算矩阵值时,不再通过平均值的方法计算像素点,而是通过正态分布的方法进行计算,其本质上是对掩模矩阵中的所有像素值进行加权平均,虽然高斯滤波克服了部分丢失边缘像素的缺陷,但由于高斯滤波在计算时,不考虑像素值的不同,因此无法完全避免。

数学公式如下:
在这里插入图片描述
API函数:

void cv::GaussianBlur ( InputArray  	src,
						OutputArray  	dst,
						Size  	ksize,
						double  	sigmaX,
						double  	sigmaY = 0,
						int  	borderType = BORDER_DEFAULT 
	) 	

作用:高斯滤波,图像高斯模糊处理

输入参数 参数定义
src 输入图像
dst 输出图像
ksize 模糊核大小,但它们必须是正的和奇数的
sigmaX X方向上的高斯核标准差
sigmaY Y方向上的高斯核标准差
borderType 边框模式用于推断图像外部的像素

中值滤波

高斯滤波对图像处理的过程与均值滤波类似,在计算矩阵值时,不再通过平均值的方法计算像素点,而是通过取中值的方法进行计算。该方法对于椒盐噪声(极大值、极小值)有很好的抑制作用,但对边缘化像素处理仍然不够完善。
如果一个3X3矩阵中为2、6、4、1、7、3、1、2、2,经过从小到大的排列1、1、2、2、2、3、4、6、7,则中值则取2,像素点赋值为2.

API函数:

void cv::medianBlur (   InputArray  	src,
						OutputArray  	dst,
						int  	ksize 
	) 	

作用:中值滤波,图像中值模糊处理

输入参数 参数定义
src 输入图像
dst 输出图像
ksize 孔径线性尺寸;它必须是奇数并且大于1

双边滤波

高斯双边滤波为了完善高斯滤波的缺点,是边缘保留的滤波方法,避免了边缘信息丢失,保留了图像轮廓的不变。相对于只考虑空间核的高斯滤波,双边滤波加入了值域核,从而将像素值相差过大的值进行保留处理,其多数用于美艳照片上。
在这里插入图片描述
API函数:

void cv::bilateralFilter (  InputArray  	src,
							OutputArray  	dst,
							int  	d,
							double  	sigmaColor,
							double  	sigmaSpace,
							int  	borderType = BORDER_DEFAULT 
	) 		

作用:双边滤波,图像双边模糊处理

输入参数 参数定义
src 输入图像
dst 输出图像
d 用于滤波的每个像素邻域的直径,直径内的所有像素均会被纳入计算。如果它是非正数,则从sigmaSpace计算
sigmaColor 决定多少差值内的像素会被计算
sigmaSpace 如果d的值大于0则无效,否则会用来计算d的值
borderType 边框模式用于推断图像外部的像素

示例

以均值滤波、高斯滤波、中值滤波、双边滤波采用5*5左右核对同一张图片进行模糊。

#include "mainwindow.h"
#include "ui_mainwindow.h"
#include <opencv2/opencv.hpp>
#include <QtDebug>

using namespace cv;

MainWindow::MainWindow(QWidget *parent) :
    QMainWindow(parent),
    ui(new Ui::MainWindow)
{
    ui->setupUi(this);
    Mat dstBlur;
    Mat dstGaussianBlur;
    Mat dstmedianBlur;
    Mat dstbilateralFilter;
    Mat src = imread("E:/OpenCV/OpenCVPicture/horse.png");
    if(src.empty()){
        qDebug()<<"can not load image...\n";
        return ;
    }
    namedWindow("srcImage",CV_WINDOW_AUTOSIZE);
    namedWindow("均值滤波",CV_WINDOW_AUTOSIZE);
    namedWindow("高斯滤波",CV_WINDOW_AUTOSIZE);
    namedWindow("中值滤波",CV_WINDOW_AUTOSIZE);
    namedWindow("双边滤波",CV_WINDOW_AUTOSIZE);
    //均值滤波
    blur(src,dstBlur,Size(5,5));
    //高斯滤波
    GaussianBlur(src,dstGaussianBlur,Size(5,5),11,11);
    //中值滤波
    medianBlur(src,dstmedianBlur,5);
    //双边滤波
    bilateralFilter(src,dstbilateralFilter,5,100,3);
    
    imshow("srcImage",src);
    imshow("均值滤波",dstBlur);
    imshow("高斯滤波",dstGaussianBlur);
    imshow("中值滤波",dstmedianBlur);
    imshow("双边滤波",dstbilateralFilter);
    waitKey(0);

原图:
在这里插入图片描述
均值滤波:
在这里插入图片描述
高斯滤波:
在这里插入图片描述
中值滤波:
在这里插入图片描述
双边滤波:
在这里插入图片描述

2015-12-17 16:49:24 liushaofang 阅读数 572

对图像进行高斯模糊的处理类,实际应用效果不错!!!

package util;

import android.annotation.SuppressLint;
import android.content.Context;
import android.graphics.Bitmap;
import android.os.Build.VERSION;
import android.renderscript.Allocation;
import android.renderscript.Allocation.MipmapControl;
import android.renderscript.Element;
import android.renderscript.RenderScript;
import android.renderscript.ScriptIntrinsicBlur;
import java.lang.reflect.Array;

public class Blur
{
  private static final String TAG = "Blur";

  @SuppressLint({"NewApi"})
  public static Bitmap fastblur(Context paramContext, Bitmap paramBitmap, int paramInt)
  {
    if (Build.VERSION.SDK_INT > 16)
    {
      Bitmap localBitmap2 = paramBitmap.copy(paramBitmap.getConfig(), true);
      RenderScript localRenderScript = RenderScript.create(paramContext);
      Allocation localAllocation1 = Allocation.createFromBitmap(localRenderScript, paramBitmap, Allocation.MipmapControl.MIPMAP_NONE, 1);
      Allocation localAllocation2 = Allocation.createTyped(localRenderScript, localAllocation1.getType());
      ScriptIntrinsicBlur localScriptIntrinsicBlur = ScriptIntrinsicBlur.create(localRenderScript, Element.U8_4(localRenderScript));
      localScriptIntrinsicBlur.setRadius(paramInt);
      localScriptIntrinsicBlur.setInput(localAllocation1);
      localScriptIntrinsicBlur.forEach(localAllocation2);
      localAllocation2.copyTo(localBitmap2);
      return localBitmap2;
    }
    Bitmap localBitmap1 = paramBitmap.copy(paramBitmap.getConfig(), true);
    if (paramInt < 1)
      return null;
    int i = localBitmap1.getWidth();
    int j = localBitmap1.getHeight();
    int[] arrayOfInt1 = new int[i * j];
    localBitmap1.getPixels(arrayOfInt1, 0, i, 0, 0, i, j);
    int k = i - 1;
    int m = j - 1;
    int n = i * j;
    int i1 = 1 + (paramInt + paramInt);
    int[] arrayOfInt2 = new int[n];
    int[] arrayOfInt3 = new int[n];
    int[] arrayOfInt4 = new int[n];
    int[] arrayOfInt5 = new int[Math.max(i, j)];
    int i2 = i1 + 1 >> 1;
    int i3 = i2 * i2;
    int[] arrayOfInt6 = new int[i3 * 256];
    for (int i4 = 0; ; i4++)
    {
      int i5 = i3 * 256;
      if (i4 >= i5)
        break;
      arrayOfInt6[i4] = (i4 / i3);
    }
    int i6 = 0;
    int i7 = 0;
    int[] arrayOfInt7 = { i1, 3 };
    int[][] arrayOfInt = (int[][])Array.newInstance(Integer.TYPE, arrayOfInt7);
    int i8 = paramInt + 1;
    for (int i9 = 0; i9 < j; i9++)
    {
      int i37 = 0;
      int i38 = 0;
      int i39 = 0;
      int i40 = 0;
      int i41 = 0;
      int i42 = 0;
      int i43 = 0;
      int i44 = 0;
      int i45 = 0;
      int i46 = -paramInt;
      if (i46 <= paramInt)
      {
        int i59 = arrayOfInt1[(i6 + Math.min(k, Math.max(i46, 0)))];
        int[] arrayOfInt13 = arrayOfInt[(i46 + paramInt)];
        arrayOfInt13[0] = ((0xFF0000 & i59) >> 16);
        arrayOfInt13[1] = ((0xFF00 & i59) >> 8);
        arrayOfInt13[2] = (i59 & 0xFF);
        int i60 = i8 - Math.abs(i46);
        i39 += i60 * arrayOfInt13[0];
        i38 += i60 * arrayOfInt13[1];
        i37 += i60 * arrayOfInt13[2];
        if (i46 > 0)
        {
          i45 += arrayOfInt13[0];
          i44 += arrayOfInt13[1];
          i43 += arrayOfInt13[2];
        }
        while (true)
        {
          i46++;
          break;
          i42 += arrayOfInt13[0];
          i41 += arrayOfInt13[1];
          i40 += arrayOfInt13[2];
        }
      }
      int i47 = paramInt;
      for (int i48 = 0; i48 < i; i48++)
      {
        arrayOfInt2[i6] = arrayOfInt6[i39];
        arrayOfInt3[i6] = arrayOfInt6[i38];
        arrayOfInt4[i6] = arrayOfInt6[i37];
        int i49 = i39 - i42;
        int i50 = i38 - i41;
        int i51 = i37 - i40;
        int[] arrayOfInt11 = arrayOfInt[((i1 + (i47 - paramInt)) % i1)];
        int i52 = i42 - arrayOfInt11[0];
        int i53 = i41 - arrayOfInt11[1];
        int i54 = i40 - arrayOfInt11[2];
        if (i9 == 0)
          arrayOfInt5[i48] = Math.min(1 + (i48 + paramInt), k);
        int i55 = arrayOfInt1[(i7 + arrayOfInt5[i48])];
        arrayOfInt11[0] = ((0xFF0000 & i55) >> 16);
        arrayOfInt11[1] = ((0xFF00 & i55) >> 8);
        arrayOfInt11[2] = (i55 & 0xFF);
        int i56 = i45 + arrayOfInt11[0];
        int i57 = i44 + arrayOfInt11[1];
        int i58 = i43 + arrayOfInt11[2];
        i39 = i49 + i56;
        i38 = i50 + i57;
        i37 = i51 + i58;
        i47 = (i47 + 1) % i1;
        int[] arrayOfInt12 = arrayOfInt[(i47 % i1)];
        i42 = i52 + arrayOfInt12[0];
        i41 = i53 + arrayOfInt12[1];
        i40 = i54 + arrayOfInt12[2];
        i45 = i56 - arrayOfInt12[0];
        i44 = i57 - arrayOfInt12[1];
        i43 = i58 - arrayOfInt12[2];
        i6++;
      }
      i7 += i;
    }
    for (int i10 = 0; i10 < i; i10++)
    {
      int i11 = 0;
      int i12 = 0;
      int i13 = 0;
      int i14 = 0;
      int i15 = 0;
      int i16 = 0;
      int i17 = 0;
      int i18 = 0;
      int i19 = 0;
      int i20 = i * -paramInt;
      int i21 = -paramInt;
      if (i21 <= paramInt)
      {
        int i35 = i10 + Math.max(0, i20);
        int[] arrayOfInt10 = arrayOfInt[(i21 + paramInt)];
        arrayOfInt10[0] = arrayOfInt2[i35];
        arrayOfInt10[1] = arrayOfInt3[i35];
        arrayOfInt10[2] = arrayOfInt4[i35];
        int i36 = i8 - Math.abs(i21);
        i13 += i36 * arrayOfInt2[i35];
        i12 += i36 * arrayOfInt3[i35];
        i11 += i36 * arrayOfInt4[i35];
        if (i21 > 0)
        {
          i19 += arrayOfInt10[0];
          i18 += arrayOfInt10[1];
          i17 += arrayOfInt10[2];
        }
        while (true)
        {
          if (i21 < m)
            i20 += i;
          i21++;
          break;
          i16 += arrayOfInt10[0];
          i15 += arrayOfInt10[1];
          i14 += arrayOfInt10[2];
        }
      }
      int i22 = i10;
      int i23 = paramInt;
      for (int i24 = 0; i24 < j; i24++)
      {
        arrayOfInt1[i22] = (0xFF000000 & arrayOfInt1[i22] | arrayOfInt6[i13] << 16 | arrayOfInt6[i12] << 8 | arrayOfInt6[i11]);
        int i25 = i13 - i16;
        int i26 = i12 - i15;
        int i27 = i11 - i14;
        int[] arrayOfInt8 = arrayOfInt[((i1 + (i23 - paramInt)) % i1)];
        int i28 = i16 - arrayOfInt8[0];
        int i29 = i15 - arrayOfInt8[1];
        int i30 = i14 - arrayOfInt8[2];
        if (i10 == 0)
          arrayOfInt5[i24] = (i * Math.min(i24 + i8, m));
        int i31 = i10 + arrayOfInt5[i24];
        arrayOfInt8[0] = arrayOfInt2[i31];
        arrayOfInt8[1] = arrayOfInt3[i31];
        arrayOfInt8[2] = arrayOfInt4[i31];
        int i32 = i19 + arrayOfInt8[0];
        int i33 = i18 + arrayOfInt8[1];
        int i34 = i17 + arrayOfInt8[2];
        i13 = i25 + i32;
        i12 = i26 + i33;
        i11 = i27 + i34;
        i23 = (i23 + 1) % i1;
        int[] arrayOfInt9 = arrayOfInt[i23];
        i16 = i28 + arrayOfInt9[0];
        i15 = i29 + arrayOfInt9[1];
        i14 = i30 + arrayOfInt9[2];
        i19 = i32 - arrayOfInt9[0];
        i18 = i33 - arrayOfInt9[1];
        i17 = i34 - arrayOfInt9[2];
        i22 += i;
      }
    }
    localBitmap1.setPixels(arrayOfInt1, 0, i, 0, 0, i, j);
    return localBitmap1;
  }
}
2015-04-21 12:55:30 smtctc 阅读数 4218

JAVA图像处理——高斯模糊

高斯模糊:
高斯模糊简单的说就是让图像的像素点取周围的像素点的平均,达到令图片模糊的效果。
当然,简单的取平均值是不太好的,因为一般与像素点距离近的像素点他们的颜色更接近,距离远的颜色偏离更大,所以可以采用正态分布曲线来取权重。
正态分布的曲线如下图:
这里写图片描述
正好符合中间点权重最高,距离中间点越远距离越低的效果,而且过度很平滑(个人这么认为)。
我们需要2维的正态分布函数,如下图:
这里写图片描述
sigma是方差,从方差的定义我们大概能推测方差越大,图像模糊的越厉害,事实也是这样。
接下来我们要决定应该取一个像素周围多少范围内的像素点来参与高斯权重值得计算呢?
我们取6sigma+1,因为大于6sigma+1的正态分布函数值(也就是权重)太小了,几乎可以忽略不计了。
接下来我们可以计算我们的高斯权重矩阵了。我们用一个二维数组来保存这些权重值。
例如我们只取周围距离一个像素的点来算平均值的画就如下图所示:
高斯权重二维矩阵
每个格子都是一个像素,对应二维矩阵内的就是当前位置像素的权重值,中间红格子是当前要计算的像素,他的颜色值=周围所有像素的 权重*颜色值之和。
根据上面的正态分布函数公式,我们写出了下面的方法求这个二维权重矩阵的所有值。

 private static float[][] getGaos(float sigma) {
        int size = (int) (6*sigma)+1;
        if(size%2==0) {
            size++;
        }
        float[][] res = new float[size][size];
        for(int y=0;y<size;y++)
            for(int x=0;x<size;x++){
                 float zhishu = -1*(((x-size/2)*(x-size/2))+((y-size/2)*(y-size/2)))/(2*sigma*sigma);
                 float a = (float) (1.0f/(2.0f*Math.PI*sigma*sigma));
                 res[x][y] = (float) (a*Math.pow(Math.E, zhishu));
                }
        return res;
    }

然后根据权重矩阵计算各个点的像素颜色值就可以了。
下面是效果图:
原图
模糊图
最后要注意边界的矩阵取不到全部的像素值时只要取能取到的值就可以了。
如果还不明白的话这里也有讲,我也是从这里参考来的:
http://www.ruanyifeng.com/blog/2012/11/gaussian_blur.html

2010-05-10 21:58:00 maozefa 阅读数 7844

阅读提示:

    《Delphi图像处理》系列以效率为侧重点,一般代码为PASCAL,核心代码采用BASM。

    《C++图像处理》系列以代码清晰,可读性为主,全部使用C++代码。

    尽可能保持二者内容一致,可相互对照。

    本文代码必须包括文章《Delphi图像处理 -- 数据类型及公用过程》中的ImageData.pas单元

 

    说明:图像高斯模糊处理代码修改次数最多,此次的修改虽然没有改变算法,但是处理流程做了修改,仅此就可以在原有基础上提高速度40%以上的。同时,这次采用了SSE浮点运算替代了原来一般汇编的定点数运算。为了方便比较,同时也不想毁去原有代码,将原代码略作修改后继续保留,此次修改的代码附在文章后面。

     我在文章《Delphi图像处理 -- 图像卷积》中,曾经介绍过利用通用的图像卷积过程对图像进行高斯模糊处理,其处理效果还不错,处理小型图像时感觉也还行,但是处理较大图像时的速度还是嫌慢,在我的P4 2.8G、1G内存的机器上对千万像素图像进行Q=3,R=5的高斯模糊处理,不包括图像装载和前期数据转换,耗时达8600ms以上,虽经几次修改,其处理速度始终得不到明显提高,主要原因还是采用通用卷积过程处理的问题:用R=5得到的卷积模板为11*11像素,一个像素有4个分量(32位ARGB),对每个象素必须作11*11*4=484个乘法、484个加法及4个除法,最后还得作4个分量是否超界的判断处理,想不慢都难啦!如果不是采用BASM定点数处理代码,其处理速度更是难以想象。

    我在网上多次查找图像高斯模糊的优化算法,不少算法和处理方式,包括代码优化还不如我的那个高斯模糊处理过程,使我很失望。前天查找其它资料时,在国外某个网站上发现介绍图像高斯模糊处理方法时似乎与常规的算法有所不同,但都没有详细的资料(因为不懂外语,很少上国外网站,但看些公式、伪代码还是行的), 经过反复琢磨,可以将其处理流程归纳如下:

    1、用给定的确定Q和长度size,计算一个size+1长的高斯分布权数数据weights: 

// 计算初始数据
for (i = -radius; i <= radius; i ++)
{
    x = i / Q;
    weights[i+radius] = exp(-x * x / 2)
}

// 求和
sum = 0
for (i = -radius; i <= radius; i ++)
{
    sum += weights[i+radius]
}

// 数据归一,即归一后的数据之和等于1
for (i = -radius; i <= radius; i ++)
{
    weights[i+radius] /= sum
}

    2、使用weights对原图像作垂直的模糊运算,即以像素(x, y)为中心,对(x, y - radius)和(x, y + radius)的像素点与weights对应的值相乘后求和得到新的像素,并写入到一个临时的图像上相应的点上(因为数据进行了归一处理,乘积和不必再作除法运算);

    3、使用weights对临时图像作水平的模糊运算,即以像素(x, y)为中心,对(x - radius, y)和(x + radius, y)的像素点与weights对应的相乘后求和得到新的像素,并写入到目标图像上相应的点上。

    处理过程结束。

    由于上面的处理流程只是对图像每个象素作了一个“十”字型的运算,使得对每个象素点的运算大大减少,模糊长度越大,减少的越多。如前面所说的Q=3、R=5的模糊运算只需要11*2*4=88个乘法、88个加法即可。

    我还是采用BASM按上面的流程作定点数运算,改进后的高斯模糊过程代码如下:

 

procedure CrossBlur(var Dest: TImageData; const Source: TImageData; Weights: Pointer; Size: Integer);
var
  height, srcStride: Integer;
  _weights: Pointer;
  dstOffset, srcOffset: Integer;
  reds, greens, blues: Integer;
asm
    push    esi
    push    edi
    push    ebx
    mov     _Weights, ecx
    mov     ecx, [edx].TImageData.Stride
    mov     srcStride, ecx
    call    _SetCopyRegs
    mov     height, edx
    mov     dstOffset, ebx
    push    esi
    push    edi
    push    edx
    push    ecx
    push    eax

    // blur col

    add     ecx, Size           // width = Source.Width
    dec     ecx
    mov     edi, _weights       // edi = weights
@@cyLoop:
    push    ecx
@@cxLoop:
    push    ecx
    push    esi
    push    edi
    xor     ebx, ebx
    mov     reds, ebx
    mov     greens, ebx
    mov     blues, ebx
    mov     ecx, Size
@@cblurLoop:
    movzx   eax, [esi].TARGBQuad.Blue
    movzx   edx, [esi].TARGBQuad.Green
    imul    eax, [edi]
    imul    edx, [edi]
    add     blues, eax
    add     greens, edx
    movzx   eax, [esi].TARGBQuad.Red
    movzx   edx, [esi].TARGBQuad.Alpha
    imul    eax, [edi]
    imul    edx, [edi]
    add     reds, eax
    add     ebx, edx
    add     edi, 4
    add     esi, srcStride
    loop    @@cblurLoop
    pop     edi
    pop     esi
    mov     eax, blues
    mov     edx, greens
    mov     ecx, reds
    shr     eax, 16
    shr     edx, 16
    shr     ecx, 16
    shr     ebx, 16
    mov     [esi].TARGBQuad.Blue, al
    mov     [esi].TARGBQuad.Green, dl
    mov     [esi].TARGBQuad.Red, cl
    mov     [esi].TARGBQuad.Alpha, bl
    add     esi, 4
    pop     ecx
    loop    @@cxLoop
    pop     ecx
    dec     height
    jnz     @@cyLoop

    pop     srcOffset
    pop     ecx
    pop     height
    pop     edi
    pop     esi

    // blur row

@@ryLoop:
    push    ecx
@@rxLoop:
    push    ecx
    push    esi
    push    edi
    xor     ebx, ebx
    mov     reds, ebx
    mov     greens, ebx
    mov     blues, ebx
    mov     ecx, Size
    mov     edi, _weights
@@rblurLoop:
    movzx   eax, [esi].TARGBQuad.Blue
    movzx   edx, [esi].TARGBQuad.Green
    imul    eax, [edi]
    imul    edx, [edi]
    add     blues, eax
    add     greens, edx
    movzx   eax, [esi].TARGBQuad.Red
    movzx   edx, [esi].TARGBQuad.Alpha
    imul    eax, [edi]
    imul    edx, [edi]
    add     reds, eax
    add     ebx, edx
    add     edi, 4
    add     esi, 4
    loop    @@rblurLoop
    pop     edi
    pop     esi
    mov     eax, blues
    mov     edx, greens
    mov     ecx, reds
    shr     eax, 16
    shr     edx, 16
    shr     ecx, 16
    shr     ebx, 16
    mov     [edi].TARGBQuad.Blue, al
    mov     [edi].TARGBQuad.Green, dl
    mov     [edi].TARGBQuad.Red, cl
    mov     [edi].TARGBQuad.Alpha, bl
    add     esi, 4
    add     edi, 4
    pop     ecx
    loop    @@rxLoop
    add     esi, srcOffset
    add     edi, dstOffset
    pop     ecx
    dec     height
    jnz     @@ryLoop
    pop     ebx
    pop     edi
    pop     esi
end;

procedure ImageGaussiabBlur(var Data: TImageData; Q: double; Radius: Integer);
var
  src: TImageData;
  fweights: array of Single;
  weights: array of Integer;
  i, size: Integer;
  fx: Double;
begin
  if Radius <= 0 then
  begin
    if Abs(Q) < 1.0 then Radius := 1
    else Radius := Round(Abs(Q)) + 2;
  end;
  size := Radius shl 1 + 1;
  SetLength(fweights, size);
  for i := 1 to Radius do
  begin
    fx := i / Q;
    fweights[Radius + i] := exp(-fx * fx / 2);
    fweights[Radius - i] := fweights[Radius + i];
  end;
  fweights[Radius] := 1.0;
  fx := 0.0;
  for i := 0 to size - 1 do
    fx := fx + fweights[i];
  SetLength(weights, size);
  for i := 0 to size - 1 do
    weights[i] := Round(fweights[i] / fx * 65536.0);
  SetLength(fweights, 0);
  src := _GetExpandData(Data, Radius);
  CrossBlur(Data, src, weights, size);
  FreeImageData(src);
end;


    用改进后的高斯模糊处理过程在我的机器上对千万像素图像进行Q=3,R=5的高斯模糊处理,不包括图像装载和前期数据转换,耗时为1390ms,处理速度确实得到了大幅度的提高。我是按32位ARGB颜色处理图像像素的,如果改为24位RGB颜色处理图像像素,耗时还可以减少,不过,RGB颜色没法处理PNG等32位像素格式的图像。

    不用模板卷积方式,而采用“十”字运算进行高斯模糊处理,效果如何呢?请看下面的简单例子代码及处理效果图:

    例子代码: 

procedure TForm1.Button3Click(Sender: TObject);
var
  bmp: TGpBitmap;
  g: TGpGraphics;
  data: TImageData;
begin
  bmp := TGpBitmap.Create('..\media\56-3.jpg');
  g := TGpGraphics.Create(Canvas.Handle);
  g.DrawImage(bmp, 0, 0);
  data := LockGpBitmap(bmp);
  ImageGaussiabBlur(Data, 3, 6);
  UnlockGpBitmap(bmp, data);
  g.DrawImage(bmp, data.Width, 0);
  g.Free;
  bmp.Free;
end;


    处理原图:

原图

    处理效果与Photoshop高斯模糊处理对比图:

效果对比图

    左上是Photoshop半径3.0高斯模糊效果图,右上是本文过程Q=3.0,R=6高斯模糊效果图。

    左下是Photoshop半径5.0高斯模糊效果图,右下是本文过程Q=5.0,R=9高斯模糊效果图。

    怎么样,效果还不错吧!

    遗憾的是我没能找到按照Q自动计算模糊半径的方法,所以处理过程给出了2个参数Q和Radius。

    下面是本次修改后的SSE代码,因原理和算法同上,只是在处理手法上有些不同:因为高斯模糊矩阵上下、左右都是对称的,因此以半径点位中心,将上下对称行(列处理时)或者左右对称列(行处理时)相加后再与高斯分布权数数据相乘,如此,除中心行(列)外,只须作以前的50%处理。

procedure CrossBlur(var Dest: TImageData; const Source: TImageData; Weights: Pointer; Radius: Integer);
var
  height, srcStride: Integer;
  dstOffset, srcOffset: Integer;
asm
    push      esi
    push      edi
    push      ebx
    push      ecx
    mov       ecx, [edx].TImageData.Stride
    mov       srcStride, ecx
    call      _SetCopyRegs
    mov       height, edx
    mov       srcOffset, eax
    mov       dstOffset, ebx
    pop       ebx
    pxor      xmm7, xmm7
    push      esi           // pst = Source.Scan0
    push      edi
    push      edx
    push      ecx

    // blur col

    mov       eax, srcStride
    mov       edx, eax
    shr       edx, 2        // width = Source.Width
    mov       edi, Radius
    shl       edi, 1
    imul      edi, eax
    add       edi, esi      // psb = pst + Radius * 2 * Source.Stride
@@cyLoop:
    push      edx
@@cxLoop:
    push      esi
    push      edi
    push      ebx
    mov       ecx, Radius
    pxor      xmm0, xmm0    // sum = 0
@@cblurLoop:
    movd      xmm1, [esi]   // for (i = 0; i < Radius; i ++)
    movd      xmm2, [edi]   // {
    punpcklbw xmm1, xmm7
    punpcklbw xmm2, xmm7
    paddw     xmm1, xmm2    //   ps = pst + psb
    punpcklwd xmm1, xmm7
    cvtdq2ps  xmm1, xmm1    //   pfs (flaot * 4) = ps (int * 4)
    mulps     xmm1, [ebx]   //   pfs *= Weights[i]
    addps     xmm0, xmm1    //   sum += pfs
    add       ebx, 16
    add       esi, eax      //   pst += Source.Stride
    sub       edi, eax      //   psb -= Source.Stride
    loop      @@cblurLoop   // }
    movd      xmm1, [esi]
    punpcklbw xmm1, xmm7
    punpcklwd xmm1, xmm7
    cvtdq2ps  xmm1, xmm1    // pfs (flaot * 4) = pst (int * 4)
    mulps     xmm1, [ebx]   // pfs *= Weights[Radius]
    addps     xmm0, xmm1    // sum += pfs
    pop       ebx
    pop       edi
    pop       esi
    cvtps2dq  xmm0, xmm0    // ps (int * 4) = sum (flaot * 4)
    packssdw  xmm0, xmm7
    packuswb  xmm0, xmm7
    movd      [esi], xmm0   // pst (byte * 4) = ps (int * 4) pask
    add       esi, 4
    add       edi, 4
    dec       edx
    jnz       @@cxLoop
    pop       edx
    dec       height
    jnz       @@cyLoop

    pop       edx
    pop       height
    pop       edi           // pd = Dest.Scan0
    pop       esi           // psl = pst
    mov       eax, Radius
    shl       eax, 1+2
    add       eax, esi      // psr = psl + Radius * 2

    // blur row

@@ryLoop:
    push      edx           // width = Dest.Width
@@rxLoop:
    push      esi
    push      ebx
    push      eax
    mov       ecx, Radius
    pxor      xmm0, xmm0    // sum = 0
@@rblurLoop:
    movd      xmm1, [esi]   // for (i = 0; i < Radius; i ++)
    movd      xmm2, [eax]   // {
    punpcklbw xmm1, xmm7
    punpcklbw xmm2, xmm7
    paddw     xmm1, xmm2    //   ps = psl + psr
    punpcklwd xmm1, xmm7
    cvtdq2ps  xmm1, xmm1    //   pfs (flaot * 4) = ps (int * 4)
    mulps     xmm1, [ebx]   //   pfs *= Weights[i]
    addps     xmm0, xmm1    //   sum += pfs
    add       ebx, 16
    add       esi, 4        //   psl ++
    sub       eax, 4        //   psr --
    loop      @@rblurLoop   // }
    movd      xmm1, [esi]
    punpcklbw xmm1, xmm7
    punpcklwd xmm1, xmm7
    cvtdq2ps  xmm1, xmm1    // pfs (flaot * 4) = psl (int * 4)
    mulps     xmm1, [ebx]   // pfs *= Weights[Radius]
    addps     xmm0, xmm1    // sum += pfs
    cvtps2dq  xmm0, xmm0    // ps (int * 4) = sum (flaot * 4)
    packssdw  xmm0, xmm7
    packuswb  xmm0, xmm7
    movd      [edi], xmm0   // pd (byte * 4) = ps (int * 4) pask
    pop       eax
    pop       ebx
    pop       esi
    add       eax, 4
    add       esi, 4
    add       edi, 4
    dec       edx
    jnz       @@rxLoop
    add       eax, srcOffset
    add       esi, srcOffset
    add       edi, dstOffset
    pop       edx
    dec       height
    jnz       @@ryLoop
    pop       ebx
    pop       edi
    pop       esi
end;

// --> st x
// <-- st e**x = 2**(x*log2(e))
function _Expon: Extended;
asm
    fldl2e              // y = x*log2e
    fmul
    fld     st(0)       // i = round(y)
    frndint
    fsub    st(1), st   // f = y - i
    fxch    st(1)       // z = 2**f
    f2xm1
    fld1
    fadd
    fscale              // result = z * 2**i
    fstp    st(1)
end;

function GetWeights(var Buffer, Weights: Pointer; Q: Single; Radius: Integer): Integer;
const
  _fcd1: Single = 0.1;
  _fc1: Single = 1.0;
  _fc2: Single = 2.0;
  _fc250: Single = 250.0;
  _fc255: Single = 255.0;
var
  R: Integer;
  v, QQ2: double;
asm
    mov     R, ecx
    mov     ecx, eax
    fld     Q
    fabs
    fcom    _fcd1
    fstsw   ax
    sahf
    jae     @@1
    fld     _fcd1
    fstp    st(1)               // if (Q < 0.1) Q = 0.1
    jmp     @@2
@@1:
    fcom    _fc250
    fstsw   ax
    sahf
    jbe     @@2
    fld     _fc250
    fstp    st(1)               // if (Q > 250) Q = 250
@@2:
    fst     Q
    fmul    Q
    fmul    _fc2
    fstp    QQ2                 // QQ2 = 2 * Q * Q
    fwait
    mov     eax, R
    test    eax, eax
    jg      @@10
    push    eax                 // if (radius <= 0)
    fld1                        // {
    fadd    Q                   //   radius = Abs(Q) + 1
    fistp   [esp].Integer
    fwait
    pop     eax
@@testRadius:                   //   while (TRUE)
    mov     R, eax              //   {
    fldz                        //     sum = 0
@@testLoop:                     //     for (R = radius; R > 0; R ++)
    fild    R                   //     {
    fld     st(0)
    fmulp   st(1), st
    fdiv    QQ2
    fchs
    call    _Expon              //       tmp = Exp(-(R * R) / (2.0 * Q * Q));
    cmp     R, eax
    jne     @@3
    fst     v                   //       if (R == radius) v = tmp
@@3:
    faddp   st(1), st(0)        //       sum += tmp
    dec     R
    jnz     @@testLoop          //     }
    fmul    _fc2                //     sum *= 2
    fadd    _fc1                //     sum += 1
    fdivr   v
    fmul    _fc255
    fistp   R
    cmp     R, 0
    je      @@4                 //     if ((INT)(v / sum * 255 + 0.5) = 0) break
    inc     eax                 //     radius ++
    jmp     @@testRadius        //   }
@@4:
    dec     eax
    jnz     @@5
    inc     eax
@@5:
    mov     R, eax              // }
@@10:
    inc     eax
    shl     eax, 4
    add     eax, 12
    push    edx
    push    ecx
    mov     edx, eax
    mov     eax, GHND
    call    GlobalAllocPtr
    pop     ecx
    pop     edx
    test    eax, eax
    jz      @@Exit
    mov     [ecx], eax          // buffer = GlobalAllocPtr(GHND, (Radius + 1) * 16 + 12)
    add     eax, 12
    and     eax, -16
    mov     [edx], eax          // weights = ((char* )buffer + 12) & 0xfffffff0
    mov     ecx, R              // ecx = radius
    mov     edx, eax            // edx = weights
    fldz                        // for (i = radius, sum = 0; i > 0; i --)
@@clacLoop:                     // {
    fild    R
    fld     st(0)
    fmulp   st(1), st
    fdiv    QQ2
    fchs
    call    _Expon
    fstp    [edx].Double        //   weights[i] = Expon(-(i * i) / (2 * Q * Q))
    fadd    [edx].Double        //   sum += weights[i]
    add     edx, 16
    dec     R
    jnz     @@clacLoop          // }
    fmul    _fc2                // sum *= 2
    fld1
    fstp    [edx].Double        // weights[radius] = 1
    fadd    [edx].Double        // sum += weights[radius]
    push    ecx
    inc     ecx
@@divLoop:                      // for (i = 0; i <= Radius; i ++)
    fld     st(0)               //   weights[i] = Round(weights[i] / sum)
    fdivr   [eax].Double
    fst     [eax].Single
    fst     [eax+4].Single
    fst     [eax+8].Single
    fstp    [eax+12].Single
    add     eax, 16
    loop    @@divLoop
    ffree   st(0)
    fwait
    pop     eax                 // return Radius
@@Exit:
end;

procedure ImageGaussiabBlur(var Data: TImageData; Q: Single; Radius: Integer);
var
  Buffer, Weights: Pointer;
  src: TImageData;
begin
  Radius := GetWeights(Buffer, Weights, Q, Radius);
  if Radius = 0 then Exit;
  if Data.AlphaFlag then
    ArgbConvertPArgb(Data);
  src := _GetExpandData(Data, Radius);
  CrossBlur(Data, src, Weights, Radius);
  FreeImageData(src);
  GlobalFreePtr(Buffer);
  if Data.AlphaFlag then
    PArgbConvertArgb(Data);
end;

    《Delphi图像处理》系列使用GDI+单元下载地址和说明见文章《GDI+ for VCL基础 -- GDI+ 与 VCL》。

    因水平有限,错误在所难免,欢迎指正和指导。邮箱地址:maozefa@hotmail.com

    这里可访问《Delphi图像处理 -- 文章索引》。

 

浅谈图像去模糊化

阅读数 683

没有更多推荐了,返回首页