把一幅图片中颜色相近的像素的分为一组，并把该组像素的颜色替换成组均值。
搞一幅比较小的bmp图像，例如128*128以内，要不然运行时间有点长。
用下面代码把像素信息导出成txt:
代码
static void BmpToTxt(string path)
{
using (StreamWriter sw = new StreamWriter("result.txt"))
{
var img = new Bitmap(path);
sw.WriteLine("dimension:{0}x{1}", img.Width, img.Height);
for (int i = 0; i < img.Width; i++)
{
for (int j = 0; j < img.Height; j++)
{
var pix = img.GetPixel(i, j);
sw.WriteLine("{0} {1} {2}", pix.R, pix.G, pix.B);
}
}
}
}
复制输出的文件为result.data，去掉第一行（图像尺寸信息）。
然后用下面的scilab代码对该数据进行meanshift：
//------------helper functions-----------------
function [feature,cls] = readdatapoint(mat,i,fieldcount,classification) //put the ith row of mat into a list
if classification == 0 then
feature = zeros(fieldcount,1);
for j=1:fieldcount
feature(j) = mat(i,j);
end
cls = -1;//not used for clustering
else
feature = zeros(fieldcount-1,1);
for j=1:fieldcount-1
feature(j) = mat(i,j);
end
cls = mat(i,fieldcount);
end
endfunction
function data = retrievedata(filepath, strfamat,classification)
if ("string" <> typeof(filepath)) then
error("filepath must be a string.");
end
fid = mopen(filepath, "r");
if(fid == -1)
error("can not open data file:" + filepath);
end
records = mfscanf(-1, fid,strfamat);
recsize = size(records);
reccount = recsize(1);
fieldcount = recsize(2);
data = list();
for i=1:reccount
[feature,cls] = readdatapoint(records,i,fieldcount, classification);
data($+1) = struct("id",i,"feature",feature,"cls",cls,"comment",1);
end
mclose(fid);
endfunction
function setsd() //reset random seed
dt = getdate();
newseed = dt(3) * dt(9)^2 + dt(10);
grand("setsd", newseed);
endfunction
//------------core functions-----------------
//compute mean shift, require global variable datapoints,n,flen
function m=calcms(y,h)
upper = zeros(flen,1);
lower = 0;
for i=1:n
xi = datapoints(i).feature;
smallkarg = (y-xi)/h;
temp = -exp(-sum(smallkarg .* smallkarg)/2)/2;
lower = lower + temp;
upper = upper + temp*xi;
end
m = lower^-1*upper-y;
endfunction
//stop if the distance between two yi and yi+1 < threshold
function y=converge(xi,h,threshold)
y = xi;
oldy = y + 10000;
while norm(y-oldy) > threshold
oldy = y;
y = y + calcms(y,h);
end
endfunction
//require global var datapoints
function clusters = runmeanshift(h,threshold,aggregatedistance)
pointcount = size(datapoints);
clusters = list();
for i=1:pointcount
pcount = i;
y = converge(datapoints(i).feature,h,threshold);
clustercount = size(clusters);
clusterid = 0; //haven't put in a cluster
clusterdistance = aggregatedistance; //a variable help to select the best cluster
for j=1:clustercount
distance = norm(clusters(j).len^-1*clusters(j).sigmay -y);//distance to cluster center
if distance <= clusterdistance then //within cluster radius
clusterdistance = distance;
clusterid=j;
end
end
if clusterid == 0 then //put in a new cluster
clusters($+1) = struct("sigmay",0, "len", 1, "elements",list(),"cls",-1);
clusters($).sigmay = y;
clusters($).elements($+1) = struct("id",i,"y",y);
else
clusters(clusterid).sigmay = clusters(clusterid).sigmay + y;
clusters(clusterid).len = clusters(clusterid).len + 1;
clusters(clusterid).elements($+1) = struct("id",i,"y",y);
end
end
endfunction
//require global clusters and datapoints
function createbmptxt(filepath)
pointnumber = size(datapoints);
pixellist = zeros(pointnumber,3);
clustercount = size(clusters);
for i=1:clustercount
clusterlength = clusters(i).len;
clustermean = clusters(i).len^-1 * clusters(i).sigmay;
clustermean = clustermean';
for j=1:clusterlength
pixellist(clusters(i).elements(j).id,:) = clustermean;
end
end
fid = mopen(filepath, "w");
if(fid == -1)
error("can not open data file:" + filepath);
end
for i=1:pointnumber
mfprintf(fid,"%d %d %d\n",pixellist(i,:));
end
mclose(fid);
endfunction
stacksize('max');
//image process code----------------------
function testimg()
datapoints = retrievedata("D:\result.data","%d %d %d",0);
flen = size(datapoints(1).feature);
flen = flen(1);
n = size(datapoints);
clusters = runmeanshift(1,50,50);
createbmptxt("D:\shiftedimg.txt");
endfunction
testimg();
可以调整runmeanshift函数的第一个参数来改变图像效果，其越大最后生成图片的颜色数就越小。
上面的代码用C#也完全可以实现，而且用C#可以使用Brahma进行GPU加速。Scilib里矩阵操作很方便，但性能好象就有点问题。
在shiftedimg.txt的第一行把图像的尺寸信息再插回去，格式和result.txt中的一样。然后用下面的代码把此文本文件转换回bmp:
代码
static void TxtToBmp(string path)
{
using (StreamReader sr = new StreamReader(path))
{
string line = sr.ReadLine();
int ind1 = line.IndexOf(":");
int ind2 = line.IndexOf("x");
int width = int.Parse(line.Substring(ind1 + 1, ind2 - ind1 - 1));
int height = int.Parse(line.Substring(ind2 + 1));
Bitmap bmp = new Bitmap(width, height);
int linecount = 0;
while ((line = sr.ReadLine()) != null)
{
string[] parts = line.Split(' ');
int r = (int)Math.Round(double.Parse(parts[0]));
int g = (int)Math.Round(double.Parse(parts[1]));
int b = (int)Math.Round(double.Parse(parts[2]));
bmp.SetPixel(linecount / height, linecount % height, Color.FromArgb(r, g, b));
linecount++;
}
bmp.Save("result.bmp");
}
}
大功告成， 最后看一下对比：
-> 