求助。。我编的程序运行时遇到未知错误，怎么改啊，谢谢众高手们了
#include <stdio.h>#include <math.h>
#include <time.h>
#include <stdlib.h>
#define n (pow(2,31)-1)
#define PI 3.14159
void read()//读取20幅图像
{
int i,j,k;
double buf[20][1024][1024]={0},buf1[1][1024][1024]={0};
double mean[1][1024][1024]={0},squaremean[1][1024][1024]={0},D_ImageMat[1][1024][1024]={0};
double ImageMat[1024][1024]={0},DImageMat[1024][1024]={0};
FILE *fp;
for (i=0;i<20;i++) //按图片的名字依次打开
{
char filename[50];
sprintf(filename,"imag_%d",i+1);
fp=fopen(filename,"r") ;
fread(buf,1,1024*1024,fp); // ImageMat(m,:,:)=fread(handle,[1024,1024],'uint16=>single');
} //打开了20 幅图片，并且分别写入20个数组里面
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
for (k=0;k<20;k++)
{
buf1[1][i][j]=buf1[1][i][j]+buf[k][1024][1024];
}
}
} //把20个数组中对应的像素中每个值都加起来
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
ImageMat[i][j]=ImageMat[i][j]+buf1[1][i][j];
}
} //因为已经相当于一幅图，所以把三维数组变成二维数组。
/////--------------求平均数-------------ImageMat[1024][1024]用于存放像素的平均数
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
ImageMat[i][j]=ImageMat[i][j]/20;
}
}
/////--------------求标准差---------------DImageMat[1024][1024]用于存放像素的标准差
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
for (k=0;k<20;k++)
{
mean[1][i][j]=(mean[1][i][j]+buf[k][i][j])/20;
squaremean[1][i][j]=(squaremean[1][i][j]+buf[k][i][j]*buf[k][i][j])/20;
D_ImageMat[1][i][j]=sqrt(squaremean[1][i][j]-mean[1][i][j]*mean[1][i][j]);
}
}
}
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
DImageMat[i][j]=DImageMat[i][j]+D_ImageMat[1][i][j];
}
}
}
void FixedPatternNoise(double ImageMat[1024][1024],double DImageMat[1024][1024])//求固定模式噪声
{
int i,j,k;
int total=0;
double sum=0.0,mDark,vDark;
for (i=0;i<1024;i++)//----求全图的平均数----
{
for (j=0;j<1024;j++)
{
sum=sum+ImageMat[i][j];
}
}
mDark=sum/(1024*1024);
//------求全图的标准差-----
sum=0.0;
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
sum=sum+DImageMat[i][j]*DImageMat[i][j];
}
}
vDark=(sqrt(sum))/1024;
//---剔除粗大误差（或 坏点，奇点噪声），并用均值替换---
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
if(fabs(ImageMat[i][j]-mDark)>6*vDark)//有大于6倍标准差的，统计个数
total++;
}
}
//----求去除粗大误差之后的固定模式噪声的均值和标准差----
for (k=1;k<=total;k++)
{
sum=0.0;
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
if(fabs(ImageMat[i][j]-mDark)>6*vDark)
ImageMat[i][j]=mDark;//大于6倍标准差的，用均值替换
}
}
sum=sum+ImageMat[i][j];
mDark=sum/(1024*1024);//替换后的均值
sum=0.0;
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
sum=(sum+ImageMat[i][j]*ImageMat[i][j])/(1024*1024);
}
vDark=sqrt(sum-mDark*mDark);//替换后的标准差
for(i=0;i<1024;i++)//再统计有超过6倍标准差的个数
{
for(j=0;j<1024;j++)
{
if (fabs(ImageMat[i][j]-mDark)>6*vDark)
total++;
}
}
if(total=0)//当没有超过6倍标准差的时候，停止循环
break;
}
}
//---------------生成对数正态分布随机数--------------------
double AverageRandom(double min,double max)//产生(min,max)之间均匀分布的随机数
{
int MINnteger = (int)(min*10000);
int MAXnteger = (int)(max*10000);
int randInteger = rand()*rand();
int diffInteger = MAXnteger - MINnteger;
int resultInteger = randInteger % diffInteger + MINnteger;
return resultInteger/10000.0;
}
double LogNormal(double x,double miu,double sigma)
//对数正态分布概率密度函数
{
return 1.0/(x*sqrt(2*PI)*sigma) * exp(-1*(log(x)-miu)*(log(x)-miu)/(2*sigma*sigma));
}
double Random_LogNormal(double miu,double sigma,double min,double max)//产生对数正态分布随机数
{
double x;
double dScope;
double y;
do
{
x = AverageRandom(min,max);
y = LogNormal(x, miu, sigma);
dScope = AverageRandom(0, LogNormal(miu,miu,sigma));
}
while( dScope > y);
return x;
}
//---------------生成正态分布随机数--------------------
double grn1(double u,double g,double *r)
{
int i,m;
double s,w,v,t;
s=65536.0;
w=2053.0;
v=13849.0;
t=0.0;
for (i=1; i<=12; i++)
{
*r=(*r)*w+v; m=(int)(*r/s);
*r=*r-m*s; t=t+(*r)/s;
}
t=u+g*(t-6.0);
return(t);
}
//------------固定模式噪声上附着的白噪声--------------------
double StdWhiteNoiseOnFPN(double ImageMat[1024][1024])
{
int i,j,k,total=0;
double StdWhiteNoiseOnFPN1,DImageMat[1024][1024],sum=0.0;
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
sum=sum+DImageMat[i][j]*DImageMat[i][j];
}
}
StdWhiteNoiseOnFPN1=(sqrt(sum))/1024;//全图的标准差
for(k=0;k<total;k++)
{
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
{
if (DImageMat[i][j]>6*StdWhiteNoiseOnFPN1)
DImageMat[i][j]=StdWhiteNoiseOnFPN1;
total++;
}
}
sum=0.0;
for (i=0;i<1024;i++)
{
for (j=0;j<1024;j++)
sum=sum+DImageMat[i][j]*DImageMat[i][j];
}
StdWhiteNoiseOnFPN1=(sqrt(sum))/1024;
if (DImageMat[i][j]>6*StdWhiteNoiseOnFPN1)
total++;
if(total=0)
break;
}
return StdWhiteNoiseOnFPN1;
}
void Poission(double m,double b[1024*1024])//生成泊松分布随机数
{
int i;
double z[1024*1024];
z[0]=m;
for(i=1;i<1024*1024;i++)
{ z[i]=fmod((16807*z[i-1]),n);
b[i]=z[i]/n;
}
b[0]=m/n;
}
void main()
{
int i,j;
int FPNplusWN[1024][1024]={0};
double mu,sigma,*r;
double mDark=0,vDark=0,StdWhiteNoiseOnFPN2;
double ImageMat[1024][1024]={0},DImageMat[1024][1024]={0};
double SimFPN[1024][1024]={0},WhiteNoiseOnFPN[1024][1024]={0};
int RadDose=50,Bth=10,T=200,lambda[1024][1024]={0},RadDoseNoise[1024][1024];
double k1=1.25,k2=0.125;
double a[1024][1024];
double b;
double c[1024*1024];
int d[1024*1024];
int k=0;
read();
FixedPatternNoise(ImageMat,DImageMat);
mu = log((mDark*mDark)/sqrt(vDark+mDark*mDark));
sigma = sqrt(log(vDark/(mDark*mDark)+1));
for(i=0;i<=1024;i++)
{
for(j=0;j<=1024;j++)
SimFPN[i][j]=Random_LogNormal(mu,sigma,mu-3*sigma,mu+3*sigma);//调用对数正态分布生成的随机数
}
StdWhiteNoiseOnFPN2=StdWhiteNoiseOnFPN(ImageMat);
mu=0;
sigma=StdWhiteNoiseOnFPN2;
for(i=0;i<=1024;i++)
{
for(j=0;j<=1024;j++)
WhiteNoiseOnFPN[i][j]=grn1(mu,sigma*sigma,r);//调用正态分布生成的随机数
}
for(i=0;i<1024;i++)
{
for(j=0;j<1024;j++)
FPNplusWN[i][j]=SimFPN[i][j]+WhiteNoiseOnFPN[i][j];
}
/////////////
for(i=0;i<=1024;i++)
{
for(j=0;j<=1024;j++)
lambda[i][j]=(RadDose-Bth)*(k1+k2*(SimFPN[i][j]-mDark));
a[i][j]=exp(-lambda[i][j]);
}
for(i=0;i<1024*1024;i++)
{
b=1;
Poission(i*123456789000,c); //调用子函数
for(j=0;j<1024*1024;j++)
{
b=b*c[j];
if(b<a[i][j])
{
d[k++]=j; //产生poisson
break;
}
}
}
k=0;
for(i=0;i<1024;i++)
{
for(j=0;j<1024;j++)
{
RadDoseNoise[i][j]=d[k];
k++;
}
}
}
