Memento

Appearance

JavaScript图像处理 ​

获取图像(Image Acquisition) ​

利用JavaScript进行图像处理的第一步是获取图像。这可以通过读取html文档里面的图像来完成。 首先，我们需要一个图像类来存储图像数据。这个图像类支持最基本的几个图像操作：新建图像，设置、取出像素点以及双线性插值采样。

js
function RGBAImage( w, h, data )
{
    this.type = 'RGBAImage';
    this.w = w;
    this.h = h;
    this.data = new Uint8Array(w*h*4);
    data && this.data.set(data);
}
 
RGBAImage.prototype.getPixel = function(x, y) {
    var idx = (y * this.w + x) * 4;
    return new Color(
        this.data[idx+0],
        this.data[idx+1],
        this.data[idx+2],
        this.data[idx+3]
    );
}
 
// bilinear sample of the image
RGBAImage.prototype.sample = function(x, y) {
    var w = this.w, h = this.h;
    var ty = Math.floor(y);
    var dy = Math.ceil(y);
 
    var lx = Math.floor(x);
    var rx = Math.ceil(x);
 
    var fx = x - lx;
    var fy = y - ty;
 
    var c = this.getPixel(lx, ty).mul((1-fy) * (1-fx))
        .add(this.getPixel(lx, dy).mul(fy * (1-fx)))
        .add(this.getPixel(rx, ty).mul((1-fy) * fx))
        .add(this.getPixel(rx, dy).mul(fy * fx));
 
    c.clamp();
 
    return c;
};
 
RGBAImage.prototype.setPixel = function(x, y, c) {
    var idx = (y * this.w + x) * 4;
    this.data[idx] = c.r;
    this.data[idx+1] = c.g;
    this.data[idx+2] = c.b;
    this.data[idx+3] = c.a;
};

此外我们还需要一个颜色类，用于进行基本的颜色操作。

js
function Color(r, g, b, a)
{
    if( arguments.length !== 4 )
    {
        this.r = this.g = this.b = this.a = 0;
    }
    else
    {
        this.r = r;
        this.g = g;
        this.b = b;
        this.a = a;
    }
}
 
Color.RED = new Color(255, 0, 0, 255);
Color.GREEN = new Color(0, 255, 0, 255);
Color.BLUE = new Color(0, 0, 255, 255);
Color.YELLOW = new Color(255, 255, 0, 255);
Color.PURPLE = new Color(255, 0, 255, 255);
Color.CYAN = new Color(0, 255, 255, 255);
Color.WHITE = new Color(255, 255, 255, 255);
Color.BLACK = new Color(0, 0, 0, 255);
Color.GRAY = new Color(128, 128, 128, 255);
 
Color.prototype.setColor = function(that)
{
    if( that != null &&
        that.constructor === Color )
    {
        this.r = that.r;
        this.g = that.g;
        this.b = that.b;
        this.a = that.a;
        return this;
    }
    else
        return null;
};
 
Color.prototype.equal = function( that ) {
    return (this.r == that.r && this.g == that.g && this.b == that.b);
}
 
Color.prototype.add = function(that) {
    return new Color(this.r + that.r, this.g + that.g, this.b + that.b, this.a + that.a);
};
 
Color.prototype.mul = function(c)
{
    return new Color(this.r * c, this.g * c, this.b * c, this.a * c);
};
 
Color.prototype.clamp = function() {
    this.r = clamp(this.r, 0, 255);
    this.g = clamp(this.g, 0, 255);
    this.b = clamp(this.b, 0, 255);
    this.a = clamp(this.a, 0, 255);
    return this;
}
 
Color.interpolate = function(c1, c2, t)
{
    return c1.mul(t).add(c2.mul(1-t));
};

有了这两个基础类之后，我们可以用以下方法来获取一个图像：

js
/* get RGBA image data from the passed image object */
RGBAImage.fromImage = function( img, cvs ) {
    var w = img.width;
    var h = img.height;
 
    // resize the canvas for drawing
    cvs.width = w;
	cvs.height = h;
	var ctx = cvs.getContext('2d');
 
    // render the image to the canvas in order to obtain image data
    ctx.drawImage(img, 0, 0);
    var imgData = ctx.getImageData(0, 0, w, h);
    var newImage = new RGBAImage(w, h, imgData.data);
    imgData = null;
 
    // clear up the canvas
    ctx.clearRect(0, 0, w, h);
    return newImage;
};

其中img是一个Image object，cvs是一个HTML canvas object。如果想要从服务器端读取一个图片并且获取其数据，可以通过以下ImageLoader来做到：

js
var ImageLoader = function(){
    this.result = undefined;
    this.loadImage = function( imgsrc, cvs ){
        var that = this;
        // create an Image object
        img = new Image();
        img.onload = function(){
            that.result = RGBAImage.fromImage(img, cvs);
            that.result.render(cvs);
        };
        img.src = imgsrc;
    };
};

有了图像，就可以把它渲染到canvas里面啦：

js
// for html canvas
RGBAImage.prototype.toImageData = function( ctx ) {
    var imgData = ctx.createImageData(this.w, this.h);
    imgData.data.set(this.data);
    return imgData;
};
 
/* render the image to the passed canvas */
RGBAImage.prototype.render = function( cvs ) {
	canvas.width = this.w;
	canvas.height = this.h;
	context.putImageData(this.toImageData(context), 0, 0);
};

以上都完成了之后，就可以进行简单的图像处理操作了。例如，将一个彩色图片转化为灰度图：

js
// utility function
// per-pixel operation
RGBAImage.prototype.map = function( f ) {
	for(var y=0;y<this.h;y++) {
		for(var x=0;x<this.w;x++) {
			this.setPixel(x, y, f(this.getPixel(x, y)));
		}
	}
	return this;
};
 
var filters = {
	'grayscale' : function( src ) {
		if( src.type === 'RGBAImage' ) {
			return src.map(function( c ) {				
				var lev = Math.round((c.r * 299 + c.g * 587 + c.b * 114) / 1000);
				c.r = c.g = c.b = lev;
				return c;
			});
		}
		else {
			throw "Not a RGBA image!";
		}
	}
};

附上一个简单的html文档：

html
<!DOCTYPE HTML>
<html lang="en">
<head>
    <meta charset="utf-8">
    <title>ImageProcJS</title>
	<script type="text/javascript" src="color.js"></script>
	<script type="text/javascript" src="image.js"></script>
	<script type="text/javascript" src="filter.js"></script>
	<script type="text/javascript" src="imageloader.js"></script>
	<script type="text/javascript" src="utils.js"></script>
	<script>
		var context, canvas;
		var loader;
		window.onload = function() {
            		loader = new ImageLoader();
			canvas = document.getElementById('cvs');
			context = canvas.getContext('2d');
			
			document.getElementById('show_btn').onclick = function() {
                		loader.loadImage('seal.jpg', canvas);
			};
			
			document.getElementById('gs_btn').onclick = function() {
                		var I = loader.result;
				I = filters.grayscale(I);
				I.render(canvas);
			}
			
			document.getElementById('files').addEventListener('change', 
			function( e ){
				handleFileSelect(e, loader, canvas);
			}, 
			false);
		}		
		
	</script>
</head>
<body>
<button id="show_btn">Show Remote Image</button>
<button id="gs_btn">Show Grayscale</button>
Upload image:
<input type="file" id="files" name="files[]" />
<output id="list"></output>
 
<p id='imageinfo'></p>
<canvas id="cvs"></canvas>
</body>
</html>

具体的示例代码可以在github上找到。

亮度对比度调整(Brightness/Contrast) ​

基本的图像处理操作之一是颜色操作，也就是根据一定的规则修改图像中像素的值。例如说改变图像的亮度、对比度，进行色彩替换，或者是对图像进行曲线操作，复杂一点的还有色彩空间矢量化等。

修改图像的亮度和对比度实际上是对像素值进行线性变换：

js
r = r * alpha + beta
g = g * alpha + beta
b = b * alpha + beta

亮度变化比较简单，就是在上面的公式中令alpha=1, beta为输入值就可以了：

js
var filters = {
    'brightness' : function( src, val ) {
        if( src.type === 'RGBAImage' ) {
            var dc = new Color(val, val, val, 0);
            return src.map(function( c ) {
                var nc = c.add(dc);
                return nc.clamp();
            });
        }
        else {
            throw "Not a RGBA image!";
        }
    }
}

用这张海豹图像作为输入，改变亮度(+128)之后，效果如下。

相应的，改变对比度时，令alpha为输入值，beta为0即可：

js
var filters = {
    'contrast' : function( src, val ) {
    if( src.type === 'RGBAImage' ) {
        var factor = Math.max((128 + val) / 128, 0);
        return src.map(function( c0 ) {
            var c = c0.mulc(factor);
            return c.clamp();
        });
    }
    else {
            throw "Not a RGBA image!";
        }
    }
}

下图是输入图片改变了对比度(+128)之后的效果

将两个操作合在一起：

js
var filters = {
    'brightnesscontrast' : function( src, alpha, beta ) {
        if( src.type === 'RGBAImage' ) {
            var factor = Math.max((128 + val) / 128, 0);
            var dc = new Color(beta, beta, beta, 0);
            return src.map(function( c0 ) {
                var c = c0.mulc(factor).add(dc);
                return c.clamp();
            });
        }
        else {
            throw "Not a RGBA image!";
        }
    }
};

得到的效果如下图(亮度+32，对比度+64)：

直方图均衡(Histogram Equalization) ​

对图像调整亮度和对比度通常可以稍微改善图像质量，但是想要获得一个较好的效果往往需要反复尝试不同的参数组合。一个自动的调整图像亮度对比度的方法可以省下不少麻烦。直方图均衡就是这样一个算法。该算法通过统计图像中像素亮度值的分布，自动地估算出一个最优的变换——变换之后像素的亮度值“均匀”地分布在0到255之间，最大化地利用了所有有效的像素值。

直方图均衡的原理很简单，就是先统计图像像素亮度值I(x,y)的概率分布函数(probabiility distribution function) p(I), 利用这个分布函数求出相应的累积分布函数(cumulative distribution function) c(I). 直方图均衡的目标是处理后图像的概率分布函数近似为均匀分布, 即 p(I') = 常数, 对应的c(I') = I' * k, 其中k也是一个常数. 具体的推倒可以参考这篇文章。 这里给出简略的说明。

考虑一个灰度图像，图像的亮度值分布p(I)可以用以下方法来计算：

js
p(I) = (亮度为I的像素)/(图像宽度x图像高度)

那么累积分布c(I)就是：

js
c(0) = p(0)
c(I) = c(I-1) + p(I)I > 0

计算直方图和累积分布曲线的代码如下：

js
// build histogram of specified image region
function histogram(img, x1, y1, x2, y2, num_bins)
{
    if( num_bins == undefined )
        num_bins = 256;
 
    var h = img.h;
    var w = img.w;
    var hist = [];
    for(var i=0;i<num_bins;i++)
        hist[i] = 0;
 
    for(var y=y1;y<y2;y++)
    {
        for(var x=x1;x<x2;x++)
        {
            var idx = (y * w + x) * 4;
            var val = Math.floor((img.data[idx] / 255.0) * (num_bins-1));
            hist[val]++;
        }
    }
 
    return hist;
}
js
// build cdf from given pdf
function buildcdf( hist, num_bins )
{
    if( num_bins == undefined )
        num_bins = 256;
 
    var cumuhist = [];
    cumuhist[0] = hist[0];
    for(var i=1;i<num_bins;i++)
        cumuhist[i] = cumuhist[i-1] + hist[i];
 
    return cumuhist;
}

一般来说，p(I)是一个形态比较复杂的曲线，例如海豹图片：

它的直方图及对应的累积分布曲线如下图所示:

我们希望处理后的图像累积分布曲线为一直线，即对于任意的亮度值k，有

js
cdf( k ) = n*k / 255

其中 n 是总像素数。

这个目标可以通过一个变换函数 k = T(I) 来做到，这个函数的意义是在新图像中亮度为 k 的像素在原图像中的亮度为 I。这个函数满足：

js
cdf( k ) = cdf( T(I) ) = n*k / 255

因为原图像中亮度为 I 的像素在新图像中亮度变为 k，所以在原图像中，亮度值小于等于 I 的像素数应该和新图像中亮度值小于等于 k 的像素数相等，因此：

js
c(I) = cdf(T(I)) = n *k / 255

稍作变换，就可以得到：

js
k = 255 * c(I)  / n

这个就是我们需要的变换函数 T(I) 。在实际计算中，可以先将累积分布函数归一化以减少计算量。

js
var filters = {
'histogram': function( src ) {
if( src.type === 'RGBAImage' ) {
        // histogram equalization, blended with orignal image
        // amount is between 0 and 1
        var h = src.h, w = src.w;

        // convert color image to grayscale image for histogram computation
        var gimg = filters.grayscale(src);

        // build histogram (pdf)
        var hist = histogram(gimg, 0, 0, w, h);

        // compute cdf
        var cumuhist = buildcdf( hist );

        // normalize cdf
        var total = cumuhist[255];
        for(var i=0;i<256;i++)
            cumuhist[i] = Math.round(cumuhist[i] / total * 255.0);  // multiply 255.0 to reduce computation

        // equalize
        return src.map(function(c0){
            var lev = Math.round((c0.r * 299 + c0.g * 587 + c0.b * 114) / 1000);  // original level
            var cI = cumuhist[lev];    // mapped level
            var ratio = cI / lev;    // use the ratio to change the image
            return c0.mulc(ratio).clamp().round();
        });
    }
    else {
        throw "Not a RGBA image!";
    }
}

};

经过处理以后，图像的直方图变成下图所示的分布：

而累积分布曲线则呈现近似线性的形态：

最终效果:

自适应直方图均衡(Adaptive Histogram Equalization) ​

直方图均衡虽然可以显著改善图像质量，但是它是一个全局修改的方法——图像的每一个像素值都用完全相同的变换函数来修改。这样做的弊端在于图像局部的细节很可能因为受到其他区域的属性的影响而丢失。例如在下图的高亮度区域，全局直方图均衡调整之后，由于亮度过高而呈现一片亮白，很难再看清楚海豹的细节纹理。

原图：

直方图均衡后：

既然全局调整效果不太理想，一个直接的想法就是进行局部调整。将图片分成若干区域，每个区域单独进行直方图均衡操作，这样各个区域就不会受到其他区域的影响而导致细节丢失了。显然，这个方法的也不是完美的。由于各个区域采用了不同的变换函数，在相邻区域的边界处的像素值在处理后几乎不可能保持连续性——只有在每个区域各自的直方图都完全一致的情况下才有可能保持连续性，而这显然和实际情况不符。

解决这个问题的办法其实也很简单，既然是由于采用了不同的变换函数而导致的不连续性，那就对变换函数插值使得最终的变换函数在整个图像上是连续的就好了(Wiki)。插值的方法跟图像双线性插值是一样的，只要将变换函数想象成一个图像就好了。把每个小区域内的变换函数想象成一个像素点，这些变换函数共同构成了一副小图像，而原图像的像素点则相当于这个小图像内的亚像素元素(sub-pixel elements)。在进行直方图均衡操作时，原图的每个像素点的变换函数由它所在的区域以及周围三个相邻区域的变换函数插值得到。插值的时候考虑该像素点离这四个区域中心的距离并以此为权重进行线性插值。如下图所示，红色像素点的变换函数T由T1, T2, T3和T4插值得到：

js
T = T1 * (1-fx) * (1-fy) + T2 * fx * (1-fy) + T3 * (1-fx) * fy + T4 * fx * fy

因为直方图均衡的变换函数本身也是线性的,所以上面的公式可以直接改写成：

js
k = 255 / n * [ c1(I) * (1-fx) * (1-fy) + c2(I) * fx * (1-fy) + c3(I) * (1-fx) * fy + c4(I) * fx * fy ]

其中c1, c2, c3, c4分别为这四个区域的累积分布函数, I是红色像素点的亮度值。

有了上面的公式之后，实现起来就变得很直接了：先将图像划分成若干小区域，在每个区域分别算好累积分布函数，然后对于每个像素点，利用上面的插值公式求出要映射到的亮度值即可。

js
var filters = {
'ahe' : function( src ) {
// find a good window size
var h = src.h, w = src.w;
    // tile size
    var tilesize = [64, 64];

    // number of bins
    var num_bins = 256;

    // number of tiles in x and y direction
    var xtiles = Math.ceil(w / tilesize[0]);
    var ytiles = Math.ceil(h / tilesize[1]);

    var cdfs = new Array(ytiles);
    for(var i=0;i<ytiles;i++)
        cdfs[i] = new Array(xtiles);

    var inv_tile_size = [1.0 / tilesize[0], 1.0 / tilesize[1]];

    var binWidth = 256 / num_bins;

    var gimg = filters.grayscale(src);

    // create histograms
    for(var i=0;i<ytiles;i++)
    {
        var y0 = i * tilesize[1];
        var y1 = Math.min(y0+tilesize[1], h);
        for(var j=0;j<xtiles;j++)
        {
            var x0 = j * tilesize[0];
            var x1 = Math.min(x0+tilesize[0], w);
            var hist = histogram(gimg, x0, y0, x1, y1, num_bins);

            var cdf = buildcdf( hist );

            var total = cdf[255];
            for(var k=0;k<256;k++)
                cdf[k] = Math.round(cdf[k] / total * 255.0);

            cdfs[i][j] = cdf;
        }
    }

    var dst = new RGBAImage(w, h);

    for(var y=0;y<h;y++)
    {
        for(var x=0;x<w;x++)
        {
            // intensity of current pixel
            var I = gimg.getPixel(x, y).r;

            // bin index
            var bin = Math.floor(I / binWidth);

            // current tile
            var tx = x * inv_tile_size[0] - 0.5;
            var ty = y * inv_tile_size[1] - 0.5;

            var xl = Math.max(Math.floor(tx), 0);
            var xr = Math.min(xl+1, xtiles-1);

            var yt = Math.max(Math.floor(ty), 0);
            var yd = Math.min(yt+1, ytiles-1);

            var fx = tx - xl;
            var fy = ty - yt;

            var cdf11 = cdfs[yt][xl][bin];
            var cdf12 = cdfs[yd][xl][bin];
            var cdf21 = cdfs[yt][xr][bin];
            var cdf22 = cdfs[yd][xr][bin];

            // bilinear interpolation
            var Iout = (1 - fx) * (1 - fy) * cdf11
                + (1 - fx) *    fy  * cdf12
                +      fx  * (1 - fy) * cdf21
                +      fx  *      fy  * cdf22;

            var ratio = Iout / I;
            var c = src.getPixel(x, y).mulc(ratio).clamp();
            dst.setPixel(x, y, c);
        }
    }

    return dst;
}

};

处理完之后效果如下，相比于全局直方图均衡来说，自适应直方图均衡的很好地保留了局部的细节信息，而且还增强了图像的纹理。但是也可以看到图像中的背景噪声相应的被放大了。

这个方法处理过后的图像直方图也是近似线性的，但是因为考虑了图像的局部特性，一部分高亮度的像素点并没有映射到很高的亮度上，而且有一部分低亮度的像素点也没有映射到很低的亮度值，所以整个曲线呈现轻微的S形。

原图
直方图均衡
自适应直方图均衡

而图像的亮度分布则显著的变均匀了:

原图
直方图均衡
自适应直方图均衡

再来看看其他图片的效果:

原图
直方图均衡
自适应直方图均衡

曲线操作(Curve Manipulation) ​

直方图均衡作为一个自动的方法虽然可以在大多数情况下获得不错的效果，但是很多时候也受限于其单一的功能而无法满足多样化的图像处理需求。尤其是在图像的艺术处理方面，直方图均衡往往并不能达到期望的效果——有时候我们需要增强图像中的高光或者是明亮的背景来突出主体，有时候又需要降低图像的整体对比度而提高某些区域的对比度来增强图像的视觉冲击力。在这些情况下，我们都需要能够直接操作图像的直方图，并且在一个不断尝试的过程中来逐步获得满意的效果。

操纵直方图最常用的方法之一就是曲线工具。曲线工具提供了一个通过修改像素值映射曲线来调整直方图的方法。对于RGB图像来说，像素值映射曲线是定义在[0, 255]之间的一个取值在[0,255]之间的函数：

js
I_new = f( I ) , 0 <= I <= 255

这个函数将每一个输入像素值映射到另外一个像素值上，从而达到直接修改图像的目的。这个曲线可以是应用于修改图像的亮度，也可以修改单个颜色通道。当用于调整图像的亮度时，它相当于引入了一个非线性的函数来对图像的亮度对比度进行调整，因为对于不同亮度的像素，进行调整时所用的参数是不一样的。相比起来，之前提到过的修改图像亮度对比度的方法(JavaScript图像处理(2) - 亮度对比度调整)采用的则是线性调整的方法，因为图像中每一个像素采用的参数都是一样的。直方图均衡(JavaScript图像处理(3) - 直方图均衡，JavaScript图像处理(4) - 自适应直方图均衡)也是一种非线性调整方法，而且是自动根据图像的统计特性进行调整的方法，其所采用的映射函数来自于图像像素亮度值的分布曲线。

具体来说，在利用曲线工具调整图像时，可以先根据曲线生成一个映射函数，将这个函数保存在一个查找表(Look-up Table)里面，调整图像时直接根据查找表来进行像素值替换或者是修改就可以了。下面的代码就是一个简单的实现：

js
// lut is the look up table defined by the input curve
filters.curve = function(src, lut, channel) { 
	switch( channel )
	{
	case 'red':
		{
			return src.map(function(c0) {
				var c = new Color(lut[c0.r], c0.g, c0.b, c0.a);
				return c.round().clamp();
			});
		}
	case 'green':
		{
			return src.map(function(c0) {
				var c = new Color(c0.r, lut[c0.g], c0.b, c0.a);
				return c.round().clamp();
			});
		}
	case 'blue':
		{
			return src.map(function(c0) {
				var c = new Color(c0.r, c0.g, lut[c0.b], c0.a);
				return c.round().clamp();
			});
		}
	case 'brightness':
	default:
		{
			return src.map(function(c0) {
				var lev = Math.round((c0.r * 299 + c0.g * 587 + c0.b * 114) / 1000);
				var bias = 1e-6;			// prevent divide by zero
				var ratio = lut[lev]/(lev + bias);
				var c = c0.mulc(ratio);
				return c.round().clamp();
			});
		}
	}
};

下面是效果图。原图如下：

调整亮度后：

单独调整红色通道，注意蓝色和绿色通道的直方图没有改变：

单独调整绿色通道：

调整蓝色通道：

另一个例子：

原图:

调整后：

减色算法(Color Reduction) ​

之前在对图像进行基于直方图的处理时提到过，RGB图像每个颜色通道有256个色阶，所以RGB颜色空间最多可以有16777216(即256^3)种颜色。这是一个巨大的空间，一个显而易见的事实是对于大多数图像而言，这个颜色空间太大了——一千六百多万种颜色，几乎可以给全高清画质(1080p)的图像每个像素点分配一种颜色！大多数情况下，一个图像里面只有远远少得多的颜色，RGB颜色空间大量的颜色并没有利用上。也就是说，我们往往可以用少得多的颜色（几十到几百之间）完全描述一幅图像。

为什么要用更少的颜色描述同一幅图像？答案是其实很简单，更少的颜色意味着可以用更少的空间存储同一幅图像，在传输图像的时候也可以在不影响内容的前提下节省流量/带宽。举个简单的例子，如果我们可以用256种颜色来描述一副图像的话，那么每个像素只需要用1个字节存储；相比之下，RGB图像每个像素需要3个字节来存储。256色的图像可以节省66%的存储空间！当然，这个假设成立的前提是我们能够找到一组颜色来很好地代表原图像里的色彩，而这并不是一个简单的事情。

根据给定的图像计算出最优颜色组合的算法称为减色算法——顾名思义，就是减少描述图像所需要的颜色数量的算法。这类算法在计算机图形学发展的早起以及互联网兴起的时候曾引起研究人员的极大兴趣，因为它们可以帮助人们利用有限的硬件、网络资源获得最好的效果或最大的收益。在计算机图形学技术日趋成熟、互联网高度发达的今天，这些算法依然在不同的应用场合发挥着重要的作用。

减色算法主要有直接量化、统计量化(population based)、颜色空间分割以及聚类等方法。此外，利用人工神经网络来进行减色操作也可以获得良好的效果。

1. 直接量化 ​

直接量化是最简单的减色算法，它只是对颜色空间进行直接的重采样来减少颜色。具体来说，直接量化对每个颜色通道单独重新采样，将每个通道的色阶从256减少到某个制定的数字。这样得到一个新的小的多的颜色空间，而原图像中的每一个像素则用在新的颜色空间中的最近邻取代：

js
...
// colors:  the desired number of colors
// it is restricted to the cube of an integer here
// levs:    color levels per channel           
var levs = Math.ceil(Math.pow(colors, 1.0/3.0));
return src.map(function(c) {
    // discretize each channel
    var r = Math.round(Math.round((c.r / 255.0) * levs) / levs * 255.0);
    var g = Math.round(Math.round((c.g / 255.0) * levs) / levs * 255.0);
    var b = Math.round(Math.round((c.b / 255.0) * levs) / levs * 255.0);
    return new Color(r, g, b, c.a);
});

直接量化的效果其实已经很不错了，用256色已经可以基本上恢复出原图的绝大部分颜色了，当然如果用更多的颜色效果会好得多。从图像的直方图也可以很清楚的看到直接量化的效果。直接量化之后的直方图只包含有限的几个峰值，它们正是直接量化时的颜色采样点。

原图:

64色图像:

256色图像:

2. 统计量化 ​

直接量化是对颜色空间进行均匀采样，这样做的最大弊端在于大量的采样点落在真实像素值分布区域之外，因而对最终的图像没有任何贡献，反而浪费了大多数采样点。这个问题和之前在做直方图均衡时提到过的问题其实是一样的，即颜色空间没有被充分利用。在做直方图均衡时，我们采样的方法是调整直方图使得累积分布曲线呈线性，从而使图像像素点的亮度值尽可能均匀地分布。这里我们可以用类似的方法来提高颜色空间中采样点的利用率——利用原图的直方图来引导采样点的选取，使得每个采样点可以大致覆盖相同数量的像素点。要做到这一点，我们需要将直方图均衡算法扩展到三维颜色空间，即在三维颜色空间中建立直方图，然后根据得到的直方图来进行有效的采样。

直接在三维空间建立直方图并不是一个十分可行的方法，原因在于三维空间的直方图有一千六百多万个(256^3)单元格，同时可以预见这个直方图必定是十分稀疏的，直接对这样一个巨大的而且稀疏的直方图进行操作只能得到一个效率低下的算法。如何才能高效地处理这个直方图呢？我们可以将这个直方图分解成三个一维的直方图——每个颜色通道建立一个直方图并且单独处理。这样我们可以很自然地得到一下这个量化算法：对每个颜色通道建立直方图，然后根据这些直方图对各个颜色通道单独采样，再利用这些采样点来组合成最终的颜色表，原图中的每个像素点用颜色表中最接近的颜色替换掉。在对颜色通道采样的时候，利用和直方图均衡类似的方法，在像素值分布多的区域进行密集采样，别的区域稀疏采样。

js
... 
var hist = colorHistogram(src, 0, 0, src.w, src.h);
var rcdf = normalizecdf( buildcdf(hist[0]) );
var gcdf = normalizecdf( buildcdf(hist[1]) );
var bcdf = normalizecdf( buildcdf(hist[2]) );
 
var levels = Math.ceil(Math.pow(colors, 1.0/3.0));
 
// get sample points using CDF
var genSamples = function(cdf) {
    var pts = [];
    var step = (1.0 - cdf[0]) / levels;
 
    for(var j=0;j<=levels;j++) {
        var p = step * j + cdf[0];
        for(var i=1;i<256;i++) {
            if( cdf[i-1] <= p && cdf[i] >= p ) {
                pts.push(i);
                break;
            }
        }
    }
    return pts;
}
 
// sample points in each channel
var rPoints = genSamples(rcdf),
    gPoints = genSamples(gcdf),
    bPoints = genSamples(bcdf);
 
// assemble the samples to a color table
return src.map(function(c) {
    // find closet r sample point
    var r = findClosest(c.r, rPoints);
 
    // find closet g sample point
    var g = findClosest(c.g, gPoints);
 
    // find closet b sample point
    var b = findClosest(c.b, bPoints); return new Color(r, g, b, c.a);
});

64色效果:

256色效果:

3. 颜色空间分割(Median-Cut) ​

前面提到利用三维直方图来优化采样点的选择，但是因为直接建立三维直方图可行性不大，所以退而求其次通过三个一维直方图来近似三维直方图。颜色空间算法则是通过另外一种方法来近似三维直方图。由于图像颜色在三维颜色空间中的分布极其稀疏，只是在某些很小的区域有比较密集的分布。如果可以定位这些小区域，进而在这些小区域上分配颜色采样点的话，就可以得到一个对不同图像具有很好针对性的采样方法。定位这些小区域的方法就是下面要描述的颜色空间分割方法——Median-Cut算法。 这个算法的本质是在颜色空间建立一棵二叉树，通过不断地细化这棵树来近似得到一个颜色三维直方图，然后再根据这棵树来分配采样点。具体来说，该算法有以下几个步骤：

  1. 获取当前图像所有颜色样本
  2. 在颜色空间中选取一个方向(R,G或者B)将这些样本等分成两部分，一般选择颜色样本包围盒的长轴方向
  3. 对分出来的两部分按照同样的方法再次进行细分，直到达到所需要的精度为止；叶子节点包含至少一个颜色样本。
  4. 完成细分后，对于每个叶子节点，考虑其包含的颜色样本。取这些颜色样本的包围盒中心处颜色作为该节点的代表颜色，即我们需要的颜色采样点。所有叶子节点的代表颜色构成最终的颜色表。

具体的实现可以参考 这段示例代码。

Median-Cut算法最大的优势在于它是一个基于图像颜色样本分布的自适应方法，不论图像中颜色样本的分布如何，总是可以生成一个和颜色样本分布匹配良好的颜色表：在颜色样本分布密集的区域内采样点分布也相对密集，其他区域则分配了较少的采样点。这是因为在进行颜色空间分割的时候，总是将一个节点内的颜色样本等分成两部分。这意味着相同数目的颜色样本总是用同样数量的采样点来代表，所以颜色样本分布密集的区域，采样点的数量自然就会多，反之则相应的比较少。

Median-Cut算法自1980年由Paul Heckbert提出之后由于其出色的表现很快流行开来，是最重要、应用最广泛的减色算法之一。常用的图像处理软件如Photoshop、GIMP等都采用了这个算法或其变种。以下是用该算法得到的效果图，可以从旁边的直方图看出，这个算法所选取的采样点分布比前面两种算法都要复杂，并且是和图像的颜色分布一致的一个采样方案。很自然的，减色后的效果也是相当的好，在只用64色的情况下也能得到远比直接量化和统计量化好的多的结果，而256色图像跟原图相比误差已经很小了。

64色效果:

256色效果:

4. 颜色聚类(k-Means Clustering) ​

颜色聚类采用了和前面的算法相同的采样策略，即根据图像颜色的分布的多少来分配采样点，其不同之处在于寻找采样点的方法。颜色聚类顾名思义，即将像素按颜色的相似程度归类。得到了颜色分类之后，从每一类中选取一个代表颜色，即可得到所需的颜色采样点来组合成最终的颜色表。聚类算法(Clustering Algorithms)是机器学习领域中最基本的算法，已经研究得非常成熟了。 在这里我们采用最简单的聚类算法——k均值聚类(k-Means clustering)来展示一下通过聚类进行减色操作的方法。关于k均值聚类的具体细节可以参考这里，下面的是利用这个方法减色得到的效果，实现的细节参考这里。值得注意的是k均值聚类的结果非常依赖于初值，因此要获得好的聚类效果，需要通过一定的呃方法来选取一个好的初始类别。从下面的结果可以看到k均值聚类减色的效果非常好，64色的误差已经很小了，而256色的结果就非常接近原图了。

64色效果:

256色效果:

5. 神经网络方法(ANN) ​

神经网络方法是一个非常有趣的方法，他采取了一个和前面四个方法不同的思路来解决减色问题。前面四个方法归结起来都来源于同一个思想，即通过在颜色空间合理地选取采样点来构造颜色表，使得减色后的图像和原图尽可能地接近。神经网络方法则有点反其道而行的意味，它从一个初始的颜色表出发，通过不断修改颜色表来改善减色效果。这个颜色表通过神经元来编码，而修改颜色表的过程则是通过神经元的反馈调节来达到。神经元的反馈调节利用所给的颜色样本来进行，对于每个颜色样本，先找到编码的颜色与之最接近的神经元，然后修改该神经元编码的颜色值，使得新的颜色和这个样本的颜色更接近，即对于某个神经元的颜色(r, g, b)：

js
r = sample.r * alpha + r * (1-alpha)
g = sample.g * alpha + g * (1-alpha)
b = sample.b * alpha + b * (1-alpha)

其中alpha是松弛参数(relaxation factor)，它的值大于0小于1。上面三个式子其实就是在样本颜色和神经元颜色之间进行线性插值，并将得到的颜色值赋予神经元。很显然这样做的效果就是使得神经元的颜色更为接近输入样本的颜色。这里的松弛参数是一个很小的值，因为我们必须缓慢的调节神经元的颜色值以保证整个神经网络逐步收敛于最小误差状态（即通过神经网络产生的图像和原图之间误差最小）。如果alpha值过大的话，调整过程中神经元的颜色值会剧烈变化，几乎可以肯定最终最终不会达到收敛。这个反馈调节的过程需要反复进行许多次， 每个颜色样本对神经网络的修正随着反馈调节次数的增加越来越小，最后当颜色样本对神经网络的修正值小于一定阈值的时候停止反馈调节。具体的实现参考这里。

64色效果:

256色效果:

几种方法的效果也可以通过所得到的图像的平均量化误差(RMSE)来比较。可以看到后面三种方法可以得到明显优于直接量化和统计量化。

方法64256
直接量化19.326710.7729
统计量化17.429111.2576
Median-Cut7.25824.1104
kmeans6.72194.2924
ANN8.20664.9518

空间滤镜(Spatial Filters) ​

前面几篇介绍了在颜色空间对图像进行处理的技术，如亮度对比度调整、直方图均衡、曲线操作、颜色量化等。另一大类图像处理技术则是在图像自身的空间——二维像素矩阵空间——来进行操作。这类图像处理技术称为空间滤镜(Spatial Filters)，通常简称为滤镜(Filters)，他们是各大图像处理软件中必备的功能之一，往往还是某些软件的杀手锏功能。这是因为设计良好的滤镜往往能够带来非凡的图像处理效果，可以毫不夸张地说，各种功能各异的滤镜拥有点石成金、化腐朽为神奇的功能。

空间滤镜，简单来说就是对于图像中的每一个像素通过其邻居像素值来调整改变自身的操作。例如最简单的3x3滤镜，它可以表示为如下操作：

js
f(x, y) = sum( w(x+dx, y+dy) * p(x+dx, y+dy) ) / sum( w(x+dx, y+dy) )

其中dx, dy = {-1, 0, 1}，w是加权函数, sum( * )表示对 * 求和。p(x, y)是在(x,y)位置的像素点的值，而f(x, y)则是经过滤镜处理之后的值。对于更大的滤镜，比如说5x5，7x7，或者更大，只需要将上面的式子中dx和dy的取值范围相应地改成[-2, 2]和[-3,3]或者[-(n-1)/2, (n-1)/2]就好了。值得注意的是，一般而言我们只考虑大小为奇数的滤镜，因为偶数大小的滤镜在实现上比奇数大小的滤镜复杂但是却没有特别的优势。此外我们可能需要在滤镜的输出值加上一个偏移量来获得满意的效果，这可以通过修改上面的式子，直接在最后加上偏移量来做到。

js
f(x, y) = sum( w(x+dx, y+dy) * p(x+dx, y+dy) ) / sum( w(x+dx, y+dy) ) + bias

根据上面的定义，空间滤镜实现起来非常容易：

js
function spatialfilter(src, f) {
// src is the source image, f is the filter
// width and height of the source image
var w = src.w, h = src.h;
// filter width and height
var wf = Math.floor((f.width - 1) / 2);
var hf = Math.floor((f.height - 1) / 2);
// filter weights
var weights = f.weights;

// filter bias
var bias = f.bias;
var invfactor = 1.0 / f.factor;

return src.map(function(c0, x, y, w, h) {
    var r = 0, g = 0, b = 0;
    for(var i=-hf, fi= 0, fidx = 0;i<=hf;i++, fi++) {
        var py = clamp(i+y, 0, h-1);
        for(var j=-wf, fj=0;j<=wf;j++, fj++, fidx++) {
            var px = clamp(j+x, 0, w-1);
            var wij = weights[fidx];
            var cij = src.getPixel(px, py);
            r += cij.r * wij;
            g += cij.g * wij;
            b += cij.b * wij;
        }
    }
    r = r * invfactor + bias;
    g = g * invfactor + bias;
    b = b * invfactor + bias;

    var c = new Color(r, g, b, c0.a);
    return c.round().clamp();
});

}

上面的代码定义了通用的滤镜操作，具体的图像处理效果还得通过构造不同的滤镜来获得。例如最简单的高斯模糊滤镜定义如下：

js
function gaussianfilter(size, sigma) {
    var weights = new Float32Array(size * size);
    // create a gaussian blur filter
 
    var cx = (size-1) * 0.5;
    var cy = (size-1) * 0.5;
    var r2 = 2.0 * sigma * sigma;
 
    var wsum = 0;
 
    for(var i= 0,idx=0;i<size;i++) {
        var dy = i - cy;
        for(var j=0;j<size;j++,idx++) {
            var dx = j - cx;
            weights[idx] = Math.exp(-(dx*dx + dy*dy) / (r2));
            wsum += weights[idx];
        }
    }
 
    return {
        width: size,
        height : size,
        factor : wsum,
        bias : 0,
        weights : weights
    };
}

对下面的海豹图像应用高斯模糊得到的效果:

对高斯模糊滤镜稍加修改，即可得到锐化滤镜：

js
function sharpenfilter(size, sigma, amount) {
    var f = gaussianfilter(size, sigma);
 
    for(var i=0;i< f.weights.length;i++) {
        f.weights[i] *= -1.0;
    }
 
    var mid = Math.floor((size * size - 1) / 2);
    f.weights[mid] = f.factor + f.weights[mid] + amount;
    f.factor = amount;
    return f;
}

锐化滤镜的效果:

图像边缘提取也可以很容易地实现，例如通过拉普拉斯算子来提取边缘：

js
function laplacian() {
    return {
        width : 3,
        height : 3,
        factor : 1.0,
        bias : 0.0,
        weights : [
            -1, -1, -1,
            -1, 8, -1,
            -1, -1, -1
        ]
    };
}

效果如下:

参考 ​

Peihong's Infinite Dimensional Space