#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#define pi 3.1415926
#define Resolution 256
#define R Resolution
#define TubeNum 256
#define TN TubeNum
#define ProjNum 200
#define PN ProjNum
#define PixelSize 2.0/256
#define PS PixelSize
#define RelativeTubeSpace 1 /* Compared to PixelSize */
#define RTS RelativeTubeSpace
#define AbsoluteTubeSpace PS*RTS
#define ATS AbsoluteTubeSpace
#define IterationNum 999
#define IN IterationNum
int main ( )
{
double begin=time(NULL);
typedef struct {
double l;
int index;
}Coef;
Coef CoefArray[2*R];
FILE * sino; //Read sinogram in
sino=fopen("Sinogram","rb");
if (sino==NULL)
{
printf("Open file failure!\n");
}
double * Sinogram=new double[PN*TN];
fread(Sinogram,sizeof(double),PN*TN,sino);
if (feof(sino)==0)
{
printf("Error!Reading sinogram has not been finished.\n");
}
fclose(sino);
FILE * ei;
ei=fopen("EI","rb");
if (ei==NULL)
{
printf("Open file failure!\n");
}
double * EI=new double[PN*TN];
fread(EI,sizeof(double),PN*TN,ei);
if (feof(ei)==0)
{
printf("Error!Reading EI has not been finished.\n");
}
fclose(ei);
int i;
double max=0;
for (i=1;i<=PN*TN;i++)
{
Sinogram[i-1]=log(EI[i-1]/Sinogram[i-1]);
if (Sinogram[i-1]>max)
{
max=Sinogram[i-1];
printf("%lf ",max);
}
}
printf("\n%lf\n",max);
double * Mu=new double[R*R];
int j;
for (j=1;j<=R*R;j++)
{
Mu[j-1]=0;
}
/*FILE * mu;
mu=fopen("Mu","rb");
if (mu==NULL)
{
printf("Open file failure!\n");
}
double * Mu=new double[R*R];
fread(Mu,sizeof(double),R*R,mu);
if (feof(mu)==0)
{
printf("Error!Reading Mu has not been finished.\n");
}
fclose(mu);*/
double * V=new double[R*R];
for (j=1;j<=R*R;j++)
{
V[j-1]=0;
}
double * W=new double[PN*TN];
FILE * coefficients;
coefficients=fopen64("Coefficients","rb");
if (coefficients==NULL)
{
printf("Open file failure!\n");
}
int count;
int m;
for (i=1;i<=PN*TN;i++)
{
//W[i-1]=1;
}
for (i=1;i<=PN*TN;i++)
{
W[i-1]=0;
fread(CoefArray,sizeof(Coef),2*R,coefficients);
if (i==1)
{
printf("%f ",256*CoefArray[1].l);
}
count=CoefArray[0].index*(-1);
for (m=1;m<=count;m++)
{
W[i-1]=W[i-1]+CoefArray[m].l;
j=CoefArray[m].index;
//Mu[j-1]=Mu[j-1]+CoefArray[m].l*Sinogram[i-1];
V[j-1]=V[j-1]+CoefArray[m].l;
}
W[i-1]=1/W[i-1];
}
FILE * result;
char ResultFileName[4]="000";
result=fopen(ResultFileName,"wb");
if (result==NULL)
{
printf("Open file failure!\n");
}
int size;
size=fwrite(&Mu[0],sizeof(double),R*R,result);
if (size<R*R)
{
printf("Error in saving sinogram!");
}
fclose(result);
double * lamda=new double[IN];
int n;
for (n=1;n<=IN;n++)
{
lamda[n-1]=1;
}
double * bEstimated=new double[PN*TN];
double * WResidual=new double[PN*TN];
double * AWResidual=new double[R*R];
for (j=1;j<=R*R;j++)
{
AWResidual[j-1]=0;
}
int hundred;
int ten;
int ge;
double IterationBegin;
for (n=1;n<=IN;n++)
{
IterationBegin=time(NULL);
hundred=n/100;
ten=(n-100*hundred)/10;
ge=n-100*hundred-10*ten;
ResultFileName[0]=char(hundred+48);
ResultFileName[1]=char(ten+48);
ResultFileName[2]=char(ge+48);
rewind(coefficients);
for (i=1;i<=PN*TN;i++)
{
bEstimated[i-1]=0;
fread(CoefArray,sizeof(Coef),2*R,coefficients);
count=CoefArray[0].index*(-1);
for (m=1;m<=count;m++)
{
j=CoefArray[m].index;
bEstimated[i-1]=bEstimated[i-1]+CoefArray[m].l*Mu[j-1];
}
if (Sinogram[i-1]-bEstimated[i-1]<0)
{
//printf("%f ",Sinogram[i-1]-bEstimated[i-1]);
}
WResidual[i-1]=W[i-1]*(Sinogram[i-1]-bEstimated[i-1]);
}
rewind(coefficients);
for (i=1;i<=PN*TN;i++)
{
fread(CoefArray,sizeof(Coef),2*R,coefficients);
count=CoefArray[0].index*(-1);
for (m=1;m<=count;m++)
{
j=CoefArray[m].index;
AWResidual[j-1]=AWResidual[j-1]+CoefArray[m].l*WResidual[i-1];
}
}
for (j=1;j<=R*R;j++)
{
int j1,j2;
if (j%R==0 )
{
j1=j/R;
j2=R;
}
else
{
j1=j/R+1;
j2=j-(j1-1)*R;
}
if (AWResidual[j-1]<0)
{
printf("%d,%d,%f,%f ",j1,j2,Mu[j-1],lamda[n-1]/V[j-1]
*AWResidual[j-1]);
}
Mu[j-1]=Mu[j-1]+(lamda[n-1]/V[j-1])*AWResidual[j-1];
}
result=fopen(ResultFileName,"wb");
if (result==NULL)
{
printf("Open file failure!\n");
}
size=fwrite(&Mu[0],sizeof(double),R*R,result);
if (size<R*R)
{
printf("Error in saving sinogram!");
}
fclose(result);
printf("Eplased time in iteration %d: %fs\n",n,time(NULL)
-IterationBegin);
}
fclose(coefficients);
printf("Total eplased time: %fs\n",time(NULL)-begin);
}
没有合适的资源?快使用搜索试试~ 我知道了~
CT经典算法SART
共1个文件
txt:1个
4星 · 超过85%的资源 74 下载量 76 浏览量
2009-10-20
09:10:22
上传
评论 7
收藏 1KB RAR 举报
温馨提示
SART代数迭代算法,CT系统中的经典算法-SART
资源推荐
资源详情
资源评论
收起资源包目录
SART代数迭代算法,CT系统中经典算法-SART.rar (1个子文件)
SART.txt 6KB
共 1 条
- 1
资源评论
- 三兑空空2012-10-09恩恩,是用C语言编的,参考中
- wshanlan2013-05-24不错,可以参考
- Flying_Fish6182011-10-14谢谢……不过有点乱……编译也有错误,但总的来说还不错……
- adazhu19882011-12-22C程序,虽然用不上,用来参考不错
stevezh
- 粉丝: 0
- 资源: 1
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功