//高斯函式
double gaussian(double x,double deta){
   return exp(-0.5*(x*x)/(deta*deta));
}
//計算當前模板係數
double BilateralWindow(double *window,int width,int height,double deta_d,double deta_r){
    double *mask=(double *)malloc(sizeof(double)*width*height);
    if(mask==NULL){
        printf("bilateral window malloc wrong\n");
        exit(0);
    }
    GaussianMask(mask,width,height,deta_d);
    double detavalue=0.0;
    double center_value=window[height/2*width+width/2];
    double k=0.0;
    double result=0.0;
    for(int j=0;j<height;j++){
        for(int i=0;i<width;i++){
            detavalue=center_value-window[j*width+i];
            mask[j*width+i]*=gaussian(detavalue,deta_r);
            k+=mask[j*width+i];
        }
    }
    for(int i=0;i<width*height;i++){
        result+=mask[i]*window[i];
    }
    free(mask);
    return result/k;
}
//雙邊濾波
void BilateralFilter(IplImage *src,IplImage *dst,int width,int height,double deta_d,double deta_r){
    double *window=(double *)malloc(sizeof(double)*width*height);
    for(int j=height/2;j<src->height-height/2;j++){
        for(int i=width/2;i<src->width-width/2;i++){
            for(int m=-height/2;m<height/2+1;m++){
                for(int n=-width/2;n<width/2+1;n++)
                    window[(m+height/2)*width+n+width/2]=cvGetReal2D(src, j+m, i+n);
            }
            double value=BilateralWindow(window,width,height,deta_d,deta_r);
            cvSetReal2D(dst, j, i, value);
        }
    }
    free(window);
}