Matlab smooth函數原理

由於項目上要用到平滑一維數組數據，參考Matlab  smooth函數轉成c++代碼

 

 

 

// x,g均為數組，具體內容略
plot(x,g);hold on,plot(x,smooth(g,50),'r');
z1 = (g1-smooth(g1,50)');  figure,plot(x,z1,'.-')

 

 

藍色為平滑前，紅色為平滑后

為了要找到缺陷，即灰度值突變很大地方，可以平滑前后相減，注意這里平滑窗寬盡量選大，選擇原則是較小甚至不影響缺陷突變的地方

    平滑前后相減

 

                     

舉例

例如g=[1 2 1 5 1 1 1 3 1 1 1]

smooth(g,5)'

ans =

1.0000    1.3333    2.0000    2.0000    1.8000    2.2000    1.4000    1.4000    1.4000    1.0000    1.0000

smooth(g,10)'

ans =

    1.0000    1.3333    2.0000    1.7143    1.7778    1.7778    1.6667    1.2857    1.4000    1.0000    1.0000

以窗寬10為例(因為除的奇數，所以最大為9)  (matlab 下標以1開始,而vc從0開始)

令b=(k-1)/2;

當i小於b時    

Y[0]=g[0] =1

Y[1]=(g[0]+g[1]+g[2])/3=4/3=1.3333

Y[2]=(g[0]+g[1]+g[2]+g[3]+g[4])/5=2

Y[3]=(g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6])/7=1.7143 

Y[3]= (g[3+(-3)]+g[3+(-2)]+g[3+(-1)]+g[3+(0)]+g[3+1]+g[3+2]+g[3+3])/(2*3+1)

那么從j=-i：i    (這里k是變動的)

Y[i]=0;

   For( j=-i;j<i+1;i++)

{

      Y[i]+=g[i+j]

}

Y[i]=y[i]/(2*i+1);

 

當i大於等於時 ，同時length(g)-i>b  (這里k是固定的)

如果n為偶數，那么應該是k=n-1；然后范圍是-(k-1)/2: (k-1)/2   b=(k-1)/2

Y[4]=(g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6] +g[7]+g[8])/9=1.7778

Y[5]=(g[1]+g[2]+g[3]+g[4]+g[5]+g[6] +g[7]+g[8] +g[9])/9=1.7778

Y[6]=(g[2]+g[3]+g[4]+g[5]+g[6]+g[7]+g[8] +g[9]+g[10])/9=1.6667

Y[6]= (g[6+(-4)]+g[6+(-3)]+g[6+(-2)]+g[6+(-1)]+g[6]+g[7]+g[8] +g[9]+g[10])/9…..

Y[i]=(g[i+(-b)]+g[i+(-b+1)]+g[i+(-b+2)]+g[i+(-b+3)]+g[i+(-b+4)]+g[i+(-b+5)]+g[i+(-b+6)] +g[i+(-b+7)]+g[i+(-b+8)])/k

-b：b即

If(i>b)||(i=b)

{

   Y[i]=0;

   For( j=-b;j<b+1;b++)

{

      Y[i]+=g[i+j]

}

Y[i]=y[i]/k;}

當i大於b  同時 length(g)-i<=b   從后往前推

Y[7]=( g[4]+ g[5]+g[6] +g[7]+g[8] +g[9]+g[10])/7=1.2857

Y[8]=( g[6] +g[7]+g[8] +g[9]+g[10])/5=1.4

Y[9]=( g[8] +g[9]+g[10])/3=1

Y[10]=(g[10] =1

當i=7時，i1=10-7=3

那么Y[7]= (g[7+(-3)]+g[7+(-2)]+g[7+(-1)]+g[7+(0)]+g[7+1]+g[7+2]+g[7+3])/(2*3+1)

當i=8時，i1=2

那么Y[8]= (g[8+(-2)]+g[8+(-1)]+g[8+(0)]+g[8+1]+g[8+2])/(2*2+1)

令 i1=(n-1)-i,那么從-i1：i1

Y[i]=0;

   For( j=-i1;j<i1+1;i1++)

{

      Y[i]+=g[i+j]

}

Y[i]=y[i]/(2*i1+1);

 

C++代碼如下，nn為數組元素個數；

 

            int b=(k-1)/2;
            vector<double> Z(nn);
            if(nn>k)
            {
                for(int i=0;i<nn;i++)
                {

                    if(i<b)
                    {
                        Z[i]=0;
                        for(int j=-i;j<i+1;j++)
                        {
                            Z[i]+=D[i+j];
                        }
                        Z[i]=Z[i]/(2*i+1);
                    }
                        
                  else  if(((i>b)||(i=b))&((nn-i)>b))
                    {
                        Z[i]=0;
                        for( int j=-b;j<b+1;j++)
                        {
                            Z[i]+=D[i+j];
                        }
                        Z[i]=Z[i]/k;
                    }
                    else
                    {
                        Z[i]=0;
                        int i1=(nn-1)-i;
                        for( int j=-i1;j<i1+1;j++)
                        {
                            Z[i]+=D[i+j];
                        }
                        Z[i]=Z[i]/(2*i1+1);
                    }
                }
            }

 

 

 

 

  平滑前

 

 平滑后

 

 平滑前后相減