写的一个非局部均匀滤波，结果由一小块出现问题
按照matlab源程序用C写的一段NLM代码，运行结果由一小块出现-2147483648，经过测试发现是在507行开始到511行的第130列到150列左右那矩阵MUL有问题。但是我不知道为啥会出现这样的计算错误，计算的average=1.#INF00，sweight=1.#INF00求大神解答啊。
新建文件夹.rar
(317.54 KB)
程序代码：#define _CRT_SECURE_NO_DEPRECATE
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <windows.h>
#include <iostream>
//定义图像尺寸
#define height 512
#define width 512
//定义相似窗口和搜索窗口维度
#define f 3
#define t 9
#define sigma 25
#define degree 4
#define h degree*sigma
//定义数组
int source[height][width];
int mir_source[height+2*f][width+2*f];
int W1[2*f+1][2*f+1];
int W2[2*f+1][2*f+1];
int MUL[2*f+1][2*f+1];
int SUB1[2*f+1][2*f+1];
int SUB2[2*f+1][2*f+1];
int result[height][width];
int psnr_c[height][width];
//定义变量
double w;
double d;
double wmax;
double sweight;
double average;
double PSNR;
//读入矩阵
void read()
{
//读入矩阵文件，并将其赋给alldata
FILE* fp=fopen("baboon_噪声.txt","r");
//FILE* source=fopen("mat.txt","r");
if(source==NULL)
{
//freopen("runhistory.txt", "a+", stdout);
printf("ERROR!\n");
return;
}
////将mat读入二维数组alldata
int i,j;
for(i=0;i<height;++i)
for(j=0;j<width;++j)
{
fscanf(fp, "%d ",&source[i][j]);//这里alldata[i][j]前一定要加&，因为这里是取地址的，和之后的输出alldata[i][j]这个数不同
}
}
//镜像函数；
void mirror()//我把行列搞反了.第三个地方出错
{
int i,j;
for(i=0;i<height+2*f;i++)
for(j=0;j<width+2*f;j++)//在mir_sourcedata[][]里进行遍历
{
if(i>(f-1)&&j>(f-1)&&i<(height+f)&&j<(width+f))//顺时针判断i,j
mir_source[i][j]=source[i-f][j-f];
else if(i<f&&j<f)//1
mir_source[i][j]=source[f-i][f-j];
else if(j>f&&i<f&&j<(width+f-1))//2
mir_source[i][j]=source[f-i][j-f];
else if(j>(width+f-1)&&i<f)//3
mir_source[i][j]=source[f-i][2*width+f-j-2];
else if(j>(width+f-1)&&i>f&&i<(height+f-1))//4
mir_source[i][j]=source[i-f][2*width+f-j-2];
else if(j>(width+f-1)&&i>(height+f-1))//5
mir_source[i][j]=source[2*height+f-i-2][2*width+f-j-2];
else if(j>f&&j<(width+f-1)&&i>(height+f-1))//6
mir_source[i][j]=source[2*height+f-i-2][j-f];
else if(j<f&&i>(height+f-1))//7
mir_source[i][j]=source[2*height+f-i-2][f-j];
else if(j<f&&i>f&&i<(height+f-1))//8
mir_source[i][j]=source[i-f][f-j];
}
}
/*
备份,正确的镜像函数
void mirror()//我把行列搞反了.第三个地方出错
{
int i,j;
for(i=0;i<height+2*f;i++)
for(j=0;j<width+2*f;j++)//在mir_sourcedata[][]里进行遍历
{
if(i>(f-1)&&j>(f-1)&&i<(height+f)&&j<(width+f))//顺时针判断i,j
mir_source[i][j]=source[i-f][j-f];
else if(i<f&&j<f)//1
mir_source[i][j]=source[f-i][f-j];
else if(j>f&&i<f&&j<(width+f-1))//2
mir_source[i][j]=source[f-i][j-f];
else if(j>(width+f-1)&&i<f)//3
mir_source[i][j]=source[f-i][2*width+f-j-2];
else if(j>(width+f-1)&&i>f&&i<(height+f-1))//4
mir_source[i][j]=source[i-f][2*width+f-j-2];
else if(j>(width+f-1)&&i>(height+f-1))//5
mir_source[i][j]=source[2*height+f-i-2][2*width+f-j-2];
else if(j>f&&j<(width+f-1)&&i>(height+f-1))//6
mir_source[i][j]=source[2*height+f-i-2][j-f];
else if(j<f&&i>(height+f-1))//7
mir_source[i][j]=source[2*height+f-i-2][f-j];
else if(j<f&&i>f&&i<(height+f-1))//8
mir_source[i][j]=source[i-f][f-j];
}
}
*/
//求最大值函数
int max_value(int a,int b)
{
int L;
if(a>=b)
L=a;
else
L=b;
return L;
}
//求最小值
int min_value(int a,int b)
{
int L;
if(a<=b)
L=a;
else
L=b;
return L;
}
//将mir_source的相应数据写入W1，W2的函数,ab代表mir_source的相应块起点
void matrix_write(int W[2*f+1][2*f+1],int a,int b)
{
int i,j;
for(i=0;i<2*f+1;i++)
for(j=0;j<2*f+1;j++)
{
W[i][j]=mir_source[i+a][j+b];//!!!!!!这里是不是有点问题
}
}
//两个矩阵相减
void matrix_sub(int A[2*f+1][2*f+1],int B[2*f+1][2*f+1],int C[2*f+1][2*f+1])
{
for(int i=0;i<2*f+1;i++)
for( int j=0;j<2*f+1;j++)
C[i][j]=abs(A[i][j]-B[i][j]);//有没有可能出现负值
}
//求矩阵的和
double matrix_sum(int A[2*f+1][2*f+1])
{
double L=0;
for(int i=0;i<2*f+1;i++)
for(int j=0;j<2*f+1;j++)
L+=A[i][j];
return L;
}
//两个矩阵相乘
/*void matrix_mul(int A[2*f+1][2*f+1],int B[2*f+1][2*f+1],int C[2*f+1][2*f+1])
{
for(int i=0;i<2*f+1;i++)
for(int j=0;j<2*f+1;j++)
for(int k=0;k<2*f+1;k++)
C[i][j]+=A[i][k]*B[k][j];
}*/
void matrix_mul_point(int A[2*f+1][2*f+1],int B[2*f+1][2*f+1],int C[2*f+1][2*f+1])
{
for(int i=0;i<2*f+1;i++)
for(int j=0;j<2*f+1;j++)
C[i][j]=A[i][j]*B[i][j];
}
//sum_psnr
double psnr_sum(int matrix_data[height][width])
{
double L=0;
for(int i=0;i<height;i++)
for(int j=0;j<width;j++)
L+=matrix_data[i][j];
return L;
}
//PSNR
void PSNR_VALUE()
{
double MSE;
double sum;
for(int i=0;i<height;i++)
for(int j=0;j<width;j++)
psnr_c[i][j]=(result[i][j]-source[i][j])*(result[i][j]-source[i][j]);
sum=psnr_sum(psnr_c);
MSE=sum/(height*width);
PSNR=10*log10(255*255/MSE);
}
//输出debug时的W1，W2，MUL
/*void print_matrix_debug(int A[2*f+1][2*f+1])
{
for(int i=0;i<2*f+1;i++)
{
for(int j=0;j<2*f+1;j++)
{
printf("%d ",A[i][j]);
}
printf("\n");
}
printf("\n\n\n");
}*/
//NLM算法.将source传入mir_source后，就可以对其进行非局部均匀滤波了。步骤按照matlab程序来
void NLM()
{
int i,j,i1,j1;
for(i=0;i<height;i++)
for(j=0;j<width;j++)
{
i1=i+f;
j1=j+f;
matrix_write(W1,i1-f,j1-f);//是对的，因为txt里行太长以至于到下一行了
//对每个点的求权值开始，都average，sweight,wmax进行初始化
w=0;
wmax=0;
average=0;
sweight=0;
d=0;
int rmin=max(i1-t,f+1);
int rmax=min(i1+t,height+f);
int smin=max(j1-t,f+1);
int smax=min(j1+t,width+f);
//freopen("runhistory.txt","a+",stdout);
//printf("第%d %d个数据rmin,rmax,smin,smax分别是 %d %d %d %d\n\n",i,j,rmin,rmax,smin,smax);
/*int rmin=max_value(i1-t,f+1);
int rmax=min_value(i1+t,height+f);
int smin=max_value(j1-t,f+1);
int smax=min_value(j1+t,width+f);*/
//在搜索块里进行求权值
for(int r=rmin;r<(rmax+1);r++)//这里加了个等号，与matlab的矩阵元素表示形式不同
{
for(int s=smin;s<(smax+1);s++)
{
if(r==i1&&s==j1)
continue;
matrix_write(W2,r-f,s-f);
//W1和W2相乘
matrix_sub(W1,W2,SUB1);
matrix_sub(W1,W2,SUB2);
matrix_mul_point(SUB1,SUB2,MUL);//乘法这错了！！因为是点乘，是各个元素对应相乘
//debug用的输出三个短矩阵.找出来的原因是矩阵出错了，也就是MUL（也不一定）就是这矩阵MUL出错了，而且应该是最后一行出错
/*if(i>507&&j>130&&j<153)
{
freopen("matrix_debug.txt","a+",stdout);
print_matrix_debug(W1);
print_matrix_debug(W2);
print_matrix_debug(MUL);
}*/
d=matrix_sum(MUL);//d的值比matlab的多一个数量级
w=exp(-d/(h*h));
//freopen("runhistory1.txt","a+",stdout);
//printf("d的值是%f\n",d);
//printf("w的值是%f\n",w);
if(w>wmax)
wmax=w;
sweight+=w;
average+=w*mir_source[r][s];
}
}
//freopen("runhistory.txt","a+",stdout);//!!!!!!!!!从506,128开始，average和sweight就开始出错了
//printf("average和sweight分别是%f %f\n",average,sweight);
average+=wmax*mir_source[i1][j1];
sweight+=wmax;
if(sweight>0)
result[i][j]=(int)(average/sweight);
else
result[i][j]=source[i][j];
/*if(i>506&&j>132&&j<153)
{
freopen("runhistory.txt","a+",stdout);
printf("average=%f\n",average);
printf("sweight=%f\n",sweight);
printf("第%d %d个数=%d\n\n",i,j,result[i][j]);
}*/
//freopen("history.txt","a+",stdout);
printf("%d %d--%d：\n",i,j,result[i][j]);
}
}
//打印滤波后的结果
void print_result()
{
FILE* fp;
fp=fopen("滤波结果.txt","w");
if(!fp)
{
printf("打开输入镜像结果错误！\n");
exit(-1);
}
for(int i=0;i<height;i++)
{
for(int j=0;j<width;j++)
{
fprintf(fp,"%d ",result[i][j]);
}
fputc('\n',fp);
}
fclose(fp);
}
//打印镜像结果
void print_mir()
{
FILE* fp;
fp=fopen("镜像结果.txt","w");
if(!fp)
{
printf("打开输入镜像结果错误！\n");
exit(-1);
}
for(int i=0;i<height+2*f;i++)
{
for(int j=0;j<width+2*f;j++)
{
fprintf(fp,"%d ",mir_source[i][j]);
}
fputc('\n',fp);
}
fclose(fp);
}
//主函数
void main()
{
clock_t start,end;
start=clock();
read();
mirror();
NLM();
PSNR_VALUE();
end=clock();
print_result();
print_mir();
printf("%d 毫秒\n",end-start);
printf("PSNR是%f db\n",PSNR);
}
