图像处理之积分图应用二(快速边缘保留滤波算法)

图像处理之积分图应用二(快速边缘保留滤波算法)

一：基本原理

传统的图像边缘保留滤波算法-如高斯双边模糊、Mean-Shift模糊等计算复杂、效率比较低，虽然有各种手段优化或者快速计算方法，当时算法相对一般码农来说理解起来比较费劲，不是一个的选择，而通过积分图像实现局部均方差的边缘保留模糊算法，计算简单而且可以做到计算量跟半径无关、跟上面提到两种边缘保留滤波(EPF)算法效率高很多。首先局部均方差滤波中计算局部均值的公式如下：

当边缘很弱的时候系数K趋近于0、该点的矫正之后的像素值就接近平均值。而当边缘很强的时候系数K趋近于1、该点的模糊之后的像素值就接近等于输入像素值。上述计算中最中意的是窗口内像素的均值与方差，计算均值可以根据积分图像很容易得到，而计算方差根据一系列的数学推导可以得到如下：

就是说可以根据积分图像通过常量次数的计算得到局部的均值与方差，让这种情况下的滤波变成一个常量时间完成的操作与窗口半径大小无关。

二：算法流程

1. 根据输入的图像计算得到积分图像，参见《图像处理之积分图像算法》

2. 根据输入的半径大小计算窗口内像素均值与方差、计算得到每个像素新的像素值

3. 循环每个个像素，重复第2步计算，得到最终的局部均方差滤波图像

三：代码实现

package com.gloomyfish.ii.demo;

import java.awt.image.BufferedImage;

/**
 * fast edge preserve filter algorithm base on integral image
 * @author zhigang jia
 * @E-mail:bfnh1998@hotmail.com
 *
 */
public class FastEPFilter extends AbstractImageOptionFilter {
	// 窗口半径大小
	private int xr;
	private int yr;
	private float sigma;
	public FastEPFilter() {
	}
	
	public void setWinsize(int radius) {
		this.xr = radius;
		this.yr = radius;
	}

	public float getSigma() {
		return sigma;
	}

	public void setSigma(float sigma) {
		this.sigma = sigma;
	}

	@Override
	public BufferedImage process(BufferedImage image) {
		long time = System.currentTimeMillis();
		int width = image.getWidth();
		int height = image.getHeight();
		// get image data
		int[] pixels = new int[width * height];
		int[] outPixels = new int[width * height];
		getRGB(image, 0, 0, width, height, pixels);
		int size = (xr * 2 + 1) * (yr * 2 + 1);
		int r = 0, g = 0, b = 0;
		float sigma2 = sigma*sigma;
		
		// per-calculate integral image
		byte[] R = new byte[width*height];
		byte[] G = new byte[width*height];
		byte[] B = new byte[width*height];
		getRGB(width, height, pixels, R, G, B);
		IntIntegralImage rii = new IntIntegralImage();
		rii.setImage(R);
		rii.process(width, height);
		IntIntegralImage gii = new IntIntegralImage();
		gii.setImage(G);
		gii.process(width, height);
		IntIntegralImage bii = new IntIntegralImage();
		bii.setImage(B);
		bii.process(width, height);
		int index = 0;
		
		for (int row = yr; row < height - yr; row++) {
			for (int col = xr; col < width - xr; col++) {
				index = row * width + col;
				r = ((pixels[index] >> 16) & 0xff);
				g = (pixels[index] >> 8) & 0xff;
				b = (pixels[index] & 0xff);
				
				int sr = rii.getBlockSum(col, row, (yr * 2 + 1), (xr * 2 + 1));
				int sg = gii.getBlockSum(col, row, (yr * 2 + 1), (xr * 2 + 1));
				int sb = bii.getBlockSum(col, row, (yr * 2 + 1), (xr * 2 + 1));
				
				float vr = rii.getBlockSquareSum(col, row, (yr * 2 + 1), (xr * 2 + 1));
				float vg = gii.getBlockSquareSum(col, row, (yr * 2 + 1), (xr * 2 + 1));
				float vb = bii.getBlockSquareSum(col, row, (yr * 2 + 1), (xr * 2 + 1));
				
				// 计算均值
				float mr = sr / size;
				float mg = sg / size;
				float mb = sb / size;
				
				// 计算方差
				float dr = (vr - (sr*sr)/size)/size;
				float dg = (vg - (sg*sg)/size)/size;
				float db = (vb - (sb*sb)/size)/size;
				
				// 计算系数K
				float kr = dr / (dr+sigma2);
				float kg = dg / (dg+sigma2);
				float kb = db / (db+sigma2);
				
				// 得到滤波后的像素值
				r = (int)((1-kr)*mr + kr*r);
				g = (int)((1-kg)*mg + kg*g);
				b = (int)((1-kb)*mb + kb*b);
				
				outPixels[row * width + col] = (0xff << 24) | (clamp(r) << 16) | (clamp(g) << 8) | clamp(b);
			}
		}
		System.out.println("FastEPFilter ->> time duration : " + (System.currentTimeMillis() - time));
		BufferedImage dest = new BufferedImage(width, height, BufferedImage.TYPE_INT_ARGB);
		setRGB(dest, 0, 0, width, height, outPixels);
		return dest;
	}
	
	/** Returns the red, green and blue planes as 3 byte arrays. */
	public void getRGB(int width, int height, int[] pixels, byte[] R, byte[] G, byte[] B) {
		int c, r, g, b;
		for (int i=0; i < width*height; i++) {
			c = pixels[i];
			r = (c&0xff0000)>>16;
			g = (c&0xff00)>>8;
			b = c&0xff;
			R[i] = (byte)r;
			G[i] = (byte)g;
			B[i] = (byte)b;
		}
	}

}

四：运行效果

半径设置为5，即窗口大小为5的时候，调整参数sigma的值即可得到此效果

其实很多磨皮的算法都是基于这个算法实现的，这个才是我想说的重点，

只有耐心看到此处才可以得到正确的答案。

五：参考：

http://imagej.net/Integral_Image_Filters#Variance

http://www.activovision.com/octavi/doku.php?id=integral_images