Mean shift演算法處理映像(Scilab & C#)

把一幅圖片中顏色相近的像素的分為一組，並把該組像素的顏色替換成組均值。

搞一幅比較小的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");
            }
        }

 

 

大功告成， 最後看一下對比：

  ->