C语言版BM3D算法

#include
#include
#include
#include
#define BlockSize1 8
#define BlockSize2 8
#define BlockStep1 3
#define BlockStep2 3
#define BlockSearch1 19
#define BlockSearch2 19
#define BlockMatch1 16
#define BlockMatch2 32
#define Threshold1 2500
#define Threshold2 400
#define Lemda1_3D 2.7//注意使用时要乘上噪声标准差
#define BetaHt 2.0
#define BetaWie 2.0
#define Width 256
#define Height 256
#define NoiseStadard 25.0
#define fillen_haar 2
#define fillen_bior 10
#define M_PI 3.1415926
///////////////////////////////////
static double KDCT[8][8];
static double KIDCT[8][8];
void DCTInit（）
{
  int ij;
  double c0c1;
  c0=1/sqrt（8.0）;
  c1=c0*sqrt（2.0）;
  for（i=0;i<8;i++）
    for（j=0;j<8;j++）
      KDCT[i][j]=cos（（2*j+1）*i*M_PI/16.0）;
  //Normalization
  for（i=0;i<8;i++）
    if（i==0）
      for（j=0;j<8;j++）
        KDCT[i][j]*=c0;
    else
      for（j=0;j<8;j++）
        KDCT[i][j]*=c1;
  for（i=0;i<8;i++）
    for（j=0;j<8;j++）
      KIDCT[i][j]=KDCT[j][i];
}
void DCT1D（double * x）//routine to perform 8 point 1-D DCT
{
  int ij;
  double y[8];//temp array
  double t;
  for（i=0;i<8;i++）
  {
    t=0.0;
    for（j=0;j<8;j++）
      t+=x[j]*KDCT[i][j];
    y[i]=t;
  }
  for（i=0;i<8;i++）
    x[i]=y[i];
}
void DCT2D（double * x）//routine to perform 8*8 point 2-D DCT
{
  int ij;
  double y[8];
  //first row direction
  for（i=0;i<8;i++）
    DCT1D（x+8*i）;
  //then column direction
  for（i=0;i<8;i++）
  {
    for（j=0;j<8;j++）
      y[j]=*（x+8*j+i）;
    DCT1D（y）;
    for（j=0;j<8;j++）
      *（x+8*j+i）=y[j];
  }
}
void IDCT1D（double * x）//routine to perform 8 point 1-D Inverse DCT
{
  int ij;
  double y[8];//temp array
  double t;
  for（i=0;i<8;i++）
  {
    t=0.0;
    for（j=0;j<8;j++）
      t+=x[j]*KIDCT[i][j];
    y[i]=t;
  }
  for（i=0;i<8;i++）
    x[i]=y[i];
}
void IDCT2D（double * x）//routine to perform 8*8 point 2-D Inverse DCT
{
  int ij;
  double y[8];
  //first row direction
  for（i=0;i<8;i++）
    IDCT1D（x+8*i）;
  //then column direction
  for（i=0;i<8;i++）
  {
    for（j=0;j<8;j++）
      y[j]=*（x+8*j+i）;
    IDCT1D（y）;
    for（j=0;j<8;j++）
      *（x+8*j+i）=y[j];
  }
}
///////////////////////////////////////////////
struct Trans_2D
{
	double Trans[BlockSize1][BlockSize1];
};
struct ValueWeight
{
	double SumValue;
	double SumWeight;
};
struct Index
{
	int IndexH;
	int IndexW;
};
///////////////////////////////////////
double WSExt（double * xint nint N）
{
  int P=2*（N-1）;
  if（n<0）
    while（n<0）
      n+=P;
  else
    if（n>=P）
      while（n>=P）
        n-=P;
  if（n>=N）
    return x[P-n];
  else
    return x[n];
}
double WAExt（double * xint nint N）
{
  int P=2*N;
  if（n<0）
    while（n<0）
      n+=P;
  else
    if（n>=P）
      while（n>=P）
        n-=P;
  if（n>=N）
    if（n==P-1）
      return 0;
    else
      return -x[P-2-n];
  else
    return x[n];
}
double HSExt（double * xint nint N）
{