SG平滑算法是由Savizkg和Golag提出来的。基于最小二乘原理的多项式平滑算法,也称卷积平滑。为啥叫多项式平滑呢?且看下去。
下面使用五点平滑算法来说明平滑过程
原理很简单如图:
把光谱一段区间的等波长间隔的5个点记为X集合,多项式平滑就是利用在波长点为Xm-2,Xm-1,Xm,Xm+1,Xm+2的数据的多项式拟合值来取代Xm,,然后依次移动,直到把光谱遍历完。
代码 order必须要比滑动窗口要小,且滑动窗口不能为偶数。
void SavitskyGolaySmoothing(float *arr,int window_size, int order, int rows, int cols){
if(window_size % 2 == 0){
throw std::logic_error("only odd window size allowed");
}
if(order >= window_size){
throw std::logic_error("Order must < window_size");
}
cv::Mat A = cv::Mat::zeros(window_size, order, CV_32FC1);
cv::Mat A_T, A_INV, B;
cv::Mat kernel;
cv::Mat result = cv::Mat::zeros(window_size, 1, CV_32FC1);
int step = int((window_size - 1)/2);
for(int i = 0; i < window_size; i++){
for(int j = 0 ; j < order; j++){
float x = pow(-step + i, j);
A.at<float>(i,j) = x;
}
}
A_T = A.t();
A_INV = (A_T * A).inv();
B = A * A_INV * A_T;
B.row(step).copyTo(kernel);
float *wrap_data = new float[step*2 + cols];
for(int row =0; row < rows; row++){
//Extend start data, size is step
for(int n = 0; n < step; n++)
{
wrap_data[n] = arr[row * cols];
}
//Copy input data
for(int col =0; col < cols; col++){
wrap_data[col + step] = arr[row * cols + col];
}
//Extend end data, size is step
for(int n = 0; n < step; n++){
wrap_data[cols + step + n] = arr[row * cols + cols -1];
}
for(int m = step; m < step + cols; m++){
for(int n = -step, j = 0; n <=step; n++, j++){
result.at<float>(0, j) = wrap_data[m + n];
}
arr[row * cols + m - step] = cv::Mat(kernel * result).at<float>(0 ,0);
}
}
delete []wrap_data;
}