VC++ OPENCV CT简单重建实现及PPT

VC++ OPENCV CT简单重建实现及PPT

I=phantom（256）; 
subplot（221）
imshow（I[]）;
title（‘256*256原始图像‘）;


[NN]=size（I）;
z=2*ceil（norm（size（I）-floor（（size（I）-1）/2）-1））+3;% radon变换默认平移点数/角度


Nt=360;  % 角度采样点数
Nd=N;   % 平移数
x=pi/180; % 角度增量
d=N/Nd; % 平移步长
theta = 1:Nt;


a=zeros（N）;


[Rxp] = radon（Itheta）;

e=floor（（z-Nd）/2）+2;
R=R（e:（Nd+e-1）:）;
R1=reshape（R256360）;

[mmnn]=size（R1）;

di=lognrnd（00.15mmnn）;
R1= 10*（R1-min（R1（:）））/（ max（R1（:））-min（R1（:）））;
I0 = 1.5e5; % incident photons; decrease this for simulating “low dose“ scans
rand（‘state‘ 0） randn（‘state‘ 0）;
yi= poissrnd（I0 * di.*exp（-R1））+3*randn（size（R1））;
if any（yi（:） == 0）
  warn（‘%d of %d values are 0 in sinogram!‘ ...
       sum（yi（:）==0） length（yi（:）））;
end
R1 = log（I0 ./ max（yi0.01））; % noisy sinogram
R1=max（R10）; 
% 显示
ff=2;
uu=22000;
v=ff*exp（R1/uu）;
subplot（222）
imagesc（R1）;
title（‘256*360有噪声平行投影‘）;
colormap（gray）
colorbar
Q=reshape（R1256360）;
g=-（Nd/2-1）:（Nd/2）;
for i=1:256
    if g（i）==0
        hl（i）=1/（4*d^2）;
    else if  mod（g（i）2）==0
            hl（i）=0;
        else
            hl（i）=（-1）/（pi^2*d^2*（g（i）^2））;
        end
    end
end
k=Nd/2:（3*Nd/2-1）; 
for m=1:Nt
 pm=Q（:m）;    
 u=conv（hlpm）;
    pm=u（k）;
    Cm=（（N-1）/2）*（1-cos（（m-1）*x）-sin（（m-1）*x））;
    for i=1:N
        for j=1:N
            Xrm=Cm+（j-1）*cos（（m-1）*x）+（i-1）*sin（（m-1）*x）;
            if Xrm<1
                n=1;
                t=abs（Xrm）-floor（abs（Xrm））;
            else
                n=floor（Xrm）;
                t=Xrm-floor（Xrm）;
            end
            if n>（Nd-1）
                n=Nd-1;
            end
            p=（1-t）*pm（n）+t*pm（n+1）;
            a（N+1-ij）=a（N+1-ij）+p;
        end
    end
end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       7306  2015-11-26 12:18  CT\CT\CT\CT.cpp

     文件    9595904  2015-11-26 13:24  CT\CT\CT\CT.ncb

     文件        868  2015-11-25 21:48  CT\CT\CT\CT.sln

    ..A..H.     12288  2015-11-26 13:24  CT\CT\CT\CT.suo

     文件       4171  2015-11-21 14:31  CT\CT\CT\CT.vcproj

     文件       1409  2015-11-26 13:24  CT\CT\CT\CT.vcproj.huxing.admin.user

     文件       5312  2015-11-21 15:03  CT\CT\CT\CT.vcxproj

     文件        927  2015-11-21 15:02  CT\CT\CT\CT.vcxproj.filters

     文件        143  2015-11-21 15:02  CT\CT\CT\CT.vcxproj.user

     文件       6510  2015-11-26 12:18  CT\CT\CT\Debug\BuildLog.htm

     文件        674  2015-11-25 11:04  CT\CT\CT\Debug\cl.command.1.tlog

     文件      23054  2015-11-25 11:04  CT\CT\CT\Debug\CL.read.1.tlog

     文件        376  2015-11-25 11:04  CT\CT\CT\Debug\CL.write.1.tlog

     文件      65536  2015-11-26 12:18  CT\CT\CT\Debug\CT.exe

     文件        405  2015-11-25 19:05  CT\CT\CT\Debug\CT.exe.embed.manifest

     文件        472  2015-11-25 19:05  CT\CT\CT\Debug\CT.exe.embed.manifest.res

     文件        387  2015-11-26 12:18  CT\CT\CT\Debug\CT.exe.intermediate.manifest

     文件     553252  2015-11-26 12:18  CT\CT\CT\Debug\CT.ilk

     文件         77  2015-11-25 11:04  CT\CT\CT\Debug\CT.lastbuildstate

     文件       2426  2015-11-25 11:04  CT\CT\CT\Debug\CT.log

     文件     180860  2015-11-26 12:18  CT\CT\CT\Debug\CT.obj

     文件    1240064  2015-11-26 12:18  CT\CT\CT\Debug\CT.pdb

     文件        194  2015-11-21 15:03  CT\CT\CT\Debug\CT_manifest.rc

     文件          2  2015-11-25 11:04  CT\CT\CT\Debug\link-cvtres.read.1.tlog

     文件          2  2015-11-25 11:04  CT\CT\CT\Debug\link-cvtres.write.1.tlog

     文件          2  2015-11-25 11:04  CT\CT\CT\Debug\link.1372-cvtres.read.1.tlog

     文件          2  2015-11-25 11:04  CT\CT\CT\Debug\link.1372-cvtres.write.1.tlog

     文件          2  2015-11-25 11:04  CT\CT\CT\Debug\link.1372.read.1.tlog

     文件          2  2015-11-25 11:04  CT\CT\CT\Debug\link.1372.write.1.tlog

     文件          2  2015-11-25 11:04  CT\CT\CT\Debug\link.2328-cvtres.read.1.tlog

............此处省略170个文件信息