DIP-HW6

实验六

一、实验目的和要求

• 目的
  • 了解双边滤波的方法
• 要求
  • 自己写C语言程序，实现所有操作，不调用其他库。

二、实验内容和原理

• 内容

  • 读入图像
  • 给图像施加椒盐噪声
  • 对含有噪声的图像进行双边滤波
  • 输出图像

• 原理

  • $I_p^{bf} = \frac{1}{W_p^{bf}}\sum_{q\in S} G_{\sigma_s}(|p-q|)G_{\sigma_r}(|Ip-Iq|Iq)$

    $W_P^{bf} = \sum_{q\in S} G_{\sigma_s}(|p-q|)G_{\sigma_r}(|I_p-I_q|)$

  • $w(i,j,k,l) = e^{(-\frac{(i-k)^2+(j-l)^2}{2\sigma_d^2}-\frac{|I(i,j)-I(k,l)|^2}{2\sigma_r^2})}$

    $I_D(i,j) = \frac{\sum_{k,l} I(k,l)w(i,j,k,l)}{\sum_{k,l}w(i,j,k,l)}$

  • 根据上述公式即可计算得到新的图像

三、源代码和分析

• 源代码

  #include<stdio.h>
  #include<stdlib.h>
  #include<math.h>
  #include<string.h>
  #include<time.h>
  #define WORD unsigned short
  #define DWORD unsigned int
  #define BYTE unsigned char
  int ph,pw,psize,lsize; 
  double const sigma_r=100,sigma_d=100;
  typedef struct BMP_FILE_HEADER
  {
  	WORD bfType;    
  	DWORD bfSize;    
  	WORD bfReserved1; 
  	WORD bfReserved2; 
  	DWORD bfOffBits;   
  }BMPFILEHEADER; 
  /* 位图信息头 */
  typedef struct BMP_INFO
  {
  	DWORD biSize;            // 信息头的大小   
  	DWORD biWidth;               // 图像的宽度     
  	DWORD biHeight;              // 图像的高度     
  	WORD biPlanes;               // 图像的位面数   
  	WORD biBitCount;             // 每个像素的位数 
  	DWORD biCompression;         // 压缩类型 
  	DWORD biSizeImage;         // 图像的大小,以字节为单位 
  	DWORD biXPelsPerMeter;       // 水平分辨率 
  	DWORD biYPelsPerMeter;       // 垂直分辨率 
  	DWORD biClrUsed;             // 使用的色彩数 
  	DWORD biClrImportant;        // 重要的颜色数 
  } BMPINF;            // 40 字节	
  BMPFILEHEADER thisBmpFileHeader;  // 定义一个 BMP 文件头的结构体
  BMPINF thisBmpInfo;               // 定义一个 BMP 文件信息结构体
  void readFileHead(FILE *fp);    //fp指向的文件的文件头读入到thisBMPFileHeader中
  void writeFileHead(FILE *fp);   //thisBMPFileHeader写入到fp指向的文件的文件头中
  void BilateralFilter(BYTE *p,BYTE* pOUT,char *name);//均值滤波，结果存在pOUT中
  void OutPutImg(BYTE *pOUT,char *name);//输出图像
  void AddNoise_salt(BYTE* p, BYTE *pNoise);//增加盐噪声
  void AddNoise_pepper(BYTE* p, BYTE *pNoise);//增加椒噪声
  double getWeight(int i,int j,int tmpi,int tmpj,int I0,int I2);//计算权重
  int main()
  {
  	FILE *fp;                    // 定义一个文件指针
    if((fp=fopen("in.bmp","rb"))==NULL)
    {
      printf("Can't open the file\n");
      return 0;
    }
    printf("open the file\n");
    readFileHead(fp);
    ph=thisBmpInfo.biHeight;//高度 
  	pw=thisBmpInfo.biWidth;//宽度 
  	psize=pw*ph*3;//图片信息字节总数 
  	lsize=pw*3;//每行信息字节总数 
  	BYTE* p=(BYTE*)malloc(psize);
  	fread(p,psize,1,fp);
  	fclose(fp);
    BYTE* pOUT=(BYTE*)malloc(psize);
    BYTE* pNoise=(BYTE*)malloc(psize);
    BilateralFilter(p,pOUT,"BilateralFiler_of_origin");  
    AddNoise_salt(p,pNoise);
    BilateralFilter(pNoise,pOUT,"BilateralFiler_of_salt");
    AddNoise_pepper(p,pNoise);
    BilateralFilter(pNoise,pOUT,"BilateralFiler_of_pepper");
  	return 0;
  }
  void OutPutImg(BYTE *pOUT,char *name)
  {
      char postfix[10]=".bmp",tmp[100];
      strcpy(tmp,name);
      strcat(tmp,postfix);
      FILE *fp=fopen(tmp,"wb");
      writeFileHead(fp);
      fwrite(pOUT,psize,1,fp);
      fclose(fp);
      printf("%s",name);
      printf(" finised!\n");
  }
  void readFileHead(FILE *fp)
  {
    fread(&thisBmpFileHeader.bfType,sizeof(thisBmpFileHeader.bfType),1,fp);
  	fread(&thisBmpFileHeader.bfSize,sizeof(thisBmpFileHeader.bfSize),1,fp);
  	fread(&thisBmpFileHeader.bfReserved1,sizeof(thisBmpFileHeader.bfReserved1),1,fp);
  	fread(&thisBmpFileHeader.bfReserved2,sizeof(thisBmpFileHeader.bfReserved2),1,fp);
  	fread(&thisBmpFileHeader.bfOffBits,sizeof(thisBmpFileHeader.bfOffBits),1,fp);
    fread(&thisBmpInfo.biSize, sizeof(thisBmpInfo.biSize), 1, fp);
  	fread(&thisBmpInfo.biWidth, sizeof(thisBmpInfo.biWidth), 1, fp);
  	fread(&thisBmpInfo.biHeight, sizeof(thisBmpInfo.biHeight), 1, fp);
  	fread(&thisBmpInfo.biPlanes, sizeof(thisBmpInfo.biPlanes), 1, fp);
  	fread(&thisBmpInfo.biBitCount, sizeof(thisBmpInfo.biBitCount), 1, fp);
  	fread(&thisBmpInfo.biCompression, sizeof(thisBmpInfo.biCompression), 1, fp);
  	fread(&thisBmpInfo.biSizeImage, sizeof(thisBmpInfo.biSizeImage), 1, fp);
  	fread(&thisBmpInfo.biXPelsPerMeter, sizeof(thisBmpInfo.biXPelsPerMeter), 1, fp);
  	fread(&thisBmpInfo.biYPelsPerMeter, sizeof(thisBmpInfo.biYPelsPerMeter), 1, fp);
  	fread(&thisBmpInfo.biClrUsed, sizeof(thisBmpInfo.biClrUsed), 1, fp);
  	fread(&thisBmpInfo.biClrImportant, sizeof(thisBmpInfo.biClrImportant), 1, fp);
  }
  void writeFileHead(FILE *fp)
  {
  	//  写入位图文件头 
  	fwrite(&thisBmpFileHeader.bfType,sizeof(thisBmpFileHeader.bfType),1,fp);
  	fwrite(&thisBmpFileHeader.bfSize,sizeof(thisBmpFileHeader.bfSize),1,fp);
  	fwrite(&thisBmpFileHeader.bfReserved1,sizeof(thisBmpFileHeader.bfReserved1),1,fp);
  	fwrite(&thisBmpFileHeader.bfReserved2,sizeof(thisBmpFileHeader.bfReserved2),1,fp);
  	fwrite(&thisBmpFileHeader.bfOffBits,sizeof(thisBmpFileHeader.bfOffBits),1,fp);
  //  写入位图信息头 
  	fwrite(&thisBmpInfo.biSize, sizeof(thisBmpInfo.biSize), 1, fp);
  	fwrite(&thisBmpInfo.biWidth, sizeof(thisBmpInfo.biWidth), 1, fp);
  	fwrite(&thisBmpInfo.biHeight, sizeof(thisBmpInfo.biHeight), 1, fp);
  	fwrite(&thisBmpInfo.biPlanes, sizeof(thisBmpInfo.biPlanes), 1, fp);
  	fwrite(&thisBmpInfo.biBitCount, sizeof(thisBmpInfo.biBitCount), 1, fp);
  	fwrite(&thisBmpInfo.biCompression, sizeof(thisBmpInfo.biCompression), 1, fp);
  	fwrite(&thisBmpInfo.biSizeImage, sizeof(thisBmpInfo.biSizeImage), 1, fp);
  	fwrite(&thisBmpInfo.biXPelsPerMeter, sizeof(thisBmpInfo.biXPelsPerMeter), 1, fp);
  	fwrite(&thisBmpInfo.biYPelsPerMeter, sizeof(thisBmpInfo.biYPelsPerMeter), 1, fp);
  	fwrite(&thisBmpInfo.biClrUsed, sizeof(thisBmpInfo.biClrUsed), 1, fp);
  	fwrite(&thisBmpInfo.biClrImportant, sizeof(thisBmpInfo.biClrImportant), 1, fp);
  }
  void AddNoise_salt(BYTE* p, BYTE *pNoise)
  {
    srand(time(0));
    int *flag=(int *)malloc(ph*pw*sizeof(int));
    memset(flag,0,sizeof(int)*ph*pw);
    int num=ph*pw*0.01;
    int i,j,tmp,_;
    for(i=0;i<ph;i++)
      for(j=0;j<pw;j++)
      {
        tmp=i*lsize+j*3;
        for(_=0;_<3;_++)
          pNoise[tmp+_]=p[tmp+_];
      }
    for(i=0;i<num;i++)
    {
      int tmpi=rand()%ph,
          tmpj=rand()%pw;
      while(flag[tmpi*pw+tmpj])
      {
        tmpi=rand()%ph;
        tmpj=rand()%pw;
      }
      flag[tmpi*pw+tmpj]=1;
      for(_=0;_<3;_++)
      {
        pNoise[tmpi*lsize+tmpj*3+_]=255;
      }
    }
    OutPutImg(pNoise,"Salt_Noise");
  }
  void AddNoise_pepper(BYTE* p, BYTE *pNoise)
  {
    srand(time(0));
    int *flag=(int *)malloc(ph*pw*sizeof(int));
    memset(flag,0,sizeof(int)*ph*pw);
    int num=ph*pw*0.01;
    int i,j,tmp,_;
    for(i=0;i<ph;i++)
      for(j=0;j<pw;j++)
      {
        tmp=i*lsize+j*3;
        for(_=0;_<3;_++)
          pNoise[tmp+_]=p[tmp+_];
      }
    for(i=0;i<num;i++)
    {
      int tmpi=rand()%ph,
          tmpj=rand()%pw;
      while(flag[tmpi*pw+tmpj])
      {
        tmpi=rand()%ph;
        tmpj=rand()%pw;
      }
      flag[tmpi*pw+tmpj]=1;
      for(_=0;_<3;_++)
      {
        pNoise[tmpi*lsize+tmpj*3+_]=0;
      }
    }
    OutPutImg(pNoise,"Pepper_Noise");
  }
  void BilateralFilter(BYTE *p,BYTE* pOUT,char *name)
  {
    int i,j,tmp,_,o;
    int fx[49]={0},fy[49]={0};
    for(i=0,_=0;i<7;i++)
      for(j=0;j<7;j++,_++)
      {
        fx[_]=i-3;
        fy[_]=j-3;
      }
    for(i=0;i<ph;i++)
    {
      for(j=0;j<pw;j++)
      {
        tmp=i*lsize+j*3;
        for(_=0;_<3;_++)
        {
          int flag=0;
          double G=0,W=0;
          for(o=0;o<49;o++)
          {
            int tmpi=i+fx[o],tmpj=j+fy[o];
            if(tmpi>=0&&tmpi<ph&&tmpj>=0&&tmpj<pw)
            {
              int t=tmpi*lsize+tmpj*3;
              double w=getWeight(i,j,tmpi,tmpj,p[tmp+_],p[t+_]);
              G+=p[t+_]*w;
              W+=w;
            }
            else   flag=1;  
          }
          if(!flag)
            pOUT[tmp+_]=(G/W > 255) ? 255 : G/W;
          else
            pOUT[tmp+_]=p[tmp+_];
        }
      }
    }
    OutPutImg(pOUT,name);
  }
  double getWeight(i,j,tmpi,tmpj,I0,I2)
  {
    double d=((double)((i-tmpi)*(i-tmpi)+(j-tmpj)*(j-tmpj)))/(2*sigma_d*sigma_d);
    double r=((double)((I0-I2)*(I0-I2)))/(2*sigma_r*sigma_r);
    return exp(-1*d-r);
  }

• 分析

  • 双边滤波函数

    利用fx，fy两个数组作为坐标偏移函数

    首先通过一个循环为坐标偏移函数赋值

    然后枚举$77$方格中的像素点，利用getWeight* 函数计算得到相应的权重，并记录累加的值。

  void BilateralFilter(BYTE *p,BYTE* pOUT,char *name)
  {
    int i,j,tmp,_,o;
    int fx[49]={0},fy[49]={0};
    for(i=0,_=0;i<7;i++)
      for(j=0;j<7;j++,_++)
      {
        fx[_]=i-3;
        fy[_]=j-3;
      }
    for(i=0;i<ph;i++)
    {
      for(j=0;j<pw;j++)
      {
        tmp=i*lsize+j*3;
        for(_=0;_<3;_++)
        {
          int flag=0;
          double G=0,W=0;
          for(o=0;o<49;o++)
          {
            int tmpi=i+fx[o],tmpj=j+fy[o];
            if(tmpi>=0&&tmpi<ph&&tmpj>=0&&tmpj<pw)
            {
              int t=tmpi*lsize+tmpj*3;
              double w=getWeight(i,j,tmpi,tmpj,p[tmp+_],p[t+_]);
              G+=p[t+_]*w;
              W+=w;
            }
            else   flag=1;  
          }
          if(!flag)
            pOUT[tmp+_]=(G/W > 255) ? 255 : G/W;
          else
            pOUT[tmp+_]=p[tmp+_];
        }
      }
    }
    OutPutImg(pOUT,name);
  }

  • 权重计算

    根据公式，直接计算各个像素点的权重

    double getWeight(i,j,tmpi,tmpj,I0,I2)
    {
      double d=((double)((i-tmpi)*(i-tmpi)+(j-tmpj)*(j-tmpj)))/(2*sigma_d*sigma_d);
      double r=((double)((I0-I2)*(I0-I2)))/(2*sigma_r*sigma_r);
      return exp(-1*d-r);
    }

心得体会

• 实验结果

  对原图的双边滤波

  

  对原图施加盐噪声

  

  对噪声干扰后的图像进行双边滤波

  

• 心得

  本次实验的代码和上次实验大体相似，区别仅在于权重的计算变得复杂一些。继续沿用上次实验坐标偏移的做法，大大简化了代码的编写难度

  通过调整两个 $\sigma$ 的值，我们可以看到lena 的模糊程度会发生变化。

  通过这次实验，我对于高斯滤波、双边滤波有了更深的理解。双边滤波综合距离和像素变化两者影响的思想及方法值得我们深思。

​

​