CT图像三维重建附源码Word格式.doc
《CT图像三维重建附源码Word格式.doc》由会员分享，可在线阅读，更多相关《CT图像三维重建附源码Word格式.doc（8页珍藏版）》请在冰点文库上搜索。
P=squeeze（ph（:
:
str2num（cell2mat（numInput））））;
%将用户的输入信息转换成数字，并从ph中得到相应的片信息P
imshow（P）%展示图片P
D=250;
%将D赋值为250,是从扇束顶点到旋转中心的像素距离。
dsensor1=2;
%正实数指定扇束传感器的间距2
F1=fanbeam（P,D,'
FanSensorSpacing'
dsensor1）;
%通过P，D等计算扇束的数据值
dsensor2=1;
%正实数指定扇束传感器的间距1
F2=fanbeam（P,D,'
dsensor2）;
dsensor3=0.25%正实数指定扇束传感器的间距0.25
[F3,sensor_pos3,fan_rot_angles3]=fanbeam（P,D,...
dsensor3）;
%通过P，D等计算扇束的数据值，并得到扇束传感器的位置sensor_pos3和旋转角度fan_rot_angles3
figure,%创建窗口
imagesc（fan_rot_angles3,sensor_pos3,F3）%根据计算出的位置和角度展示F3的图片
colormap（hot）;
%设置色图为hot
colorbar;
%显示色栏
xlabel（'
FanRotationAngle（degrees）'
）%定义x坐标轴
ylabel（'
FanSensorPosition（degrees）'
）%定义y坐标轴
output_size=max（size（P））;
%得到P维数的最大值，并赋值给output_size
Ifan1=ifanbeam（F1,D,...
'
dsensor1,'
OutputSize'
output_size）;
%根据扇束投影数据F1及D等数据重建图像
figure,imshow（Ifan1）%创建窗口，并展示图片Ifan1
title（'
图一'
）;
disp（'
图一和原图的性噪比为：
result=psnr1（Ifan1,P）;
Ifan2=ifanbeam（F2,D,...
dsensor2,'
%根据扇束投影数据F2及D等数据重建图像
figure,imshow（Ifan2）%创建窗口，并展示图片Ifan2
图二和原图的性噪比为：
result=psnr1（Ifan2,P）;
图二'
Ifan3=ifanbeam（F3,D,...
dsensor3,'
%根据扇束投影数据F3及D等数据重建图像
figure,imshow（Ifan3）%创建窗口，并展示图片Ifan3
图三'
图三和原图的性噪比为：
result=psnr1（Ifan3,P）;
function[p,ellipse]=phantom3d（varargin）
%PHANTOM3DThree-dimensionalanalogueofMATLABShepp-Loganphantom
%P=PHANTOM3D（DEF,N）generatesa3Dheadphantomthatcan
%beusedtotest3-Dreconstructionalgorithms.
%
%DEFisastringthatspecifiesthetypeofheadphantomtogenerate.
%Validvaluesare:
%
%'
Shepp-Logan'
Atestimageusedwidelybyresearchersin
%tomography
ModifiedShepp-Logan'
（default）AvariantoftheShepp-Loganphantom
%inwhichthecontrastisimprovedforbetter
%visualperception.
%NisascalarthatspecifiesthegridsizeofP.
%Ifyouomittheargument,Ndefaultsto64.
%
%P=PHANTOM3D（E,N）generatesauser-definedphantom,whereeachrow
%ofthematrixEspecifiesanellipsoidintheimage.Ehastencolumns,
%witheachcolumncontainingadifferentparameterfortheellipsoids:
%
%Column1:
Atheadditiveintensityvalueoftheellipsoid
%Column2:
athelengthofthexsemi-axisoftheellipsoid
%Column3:
bthelengthoftheysemi-axisoftheellipsoid
%Column4:
cthelengthofthezsemi-axisoftheellipsoid
%Column5:
x0thex-coordinateofthecenteroftheellipsoid
%Column6:
y0they-coordinateofthecenteroftheellipsoid
%Column7:
z0thez-coordinateofthecenteroftheellipsoid
%Column8:
phiphiEulerangle（indegrees）（rotationaboutz-axis）
%Column9:
thetathetaEulerangle（indegrees）（rotationaboutx-axis）
%Column10:
psipsiEulerangle（indegrees）（rotationaboutz-axis）
%Forpurposesofgeneratingthephantom,thedomainsforthex-,y-,and
%z-axesspan[-1,1].Columns2through7mustbespecifiedinterms
%ofthisrange.
%[P,E]=PHANTOM3D（...）returnsthematrixEusedtogeneratethephantom.
%ClassSupport
%-------------
%Allinputsmustbeofclassdouble.Alloutputsareofclassdouble.
%Remarks
%-------
%Foranygivenvoxelintheoutputimage,thevoxel'
svalueisequaltothe
%sumoftheadditiveintensityvaluesofallellipsoidsthatthevoxelisa
%partof.Ifavoxelisnotpartofanyellipsoid,itsvalueis0.
%TheadditiveintensityvalueAforanellipsoidcanbepositiveornegative;
%ifitisnegative,theellipsoidwillbedarkerthanthesurroundingpixels.
%Notethat,dependingonthevaluesofA,somevoxelsmayhavevaluesoutside
%therange[0,1].
%
%Example
%ph=phantom3d（128）;
%figure,imshow（squeeze（ph（64,:
）））
%Copyright2005MatthiasChristianSchabel（matthias@stanfordalumni.org）
%UniversityofUtahDepartmentofRadiology
%UtahCenterforAdvancedImagingResearch
%729ArapeenDrive
%SaltLakeCity,UT84108-1218
%ThiscodeisreleasedundertheGnuPublicLicense（GPL）.Formoreinformation,
%see:
http:
//www.gnu.org/copyleft/gpl.html
%Portionsofthiscodearebasedonphantom.m,copyrightedbytheMathworks
[ellipse,n]=parse_inputs（varargin{:
}）;
p=zeros（[nnn]）;
rng=（（0:
n-1）-（n-1）/2）/（（n-1）/2）;
[x,y,z]=meshgrid（rng,rng,rng）;
coord=[flatten（x）;
flatten（y）;
flatten（z）];
p=flatten（p）;
fork=1:
size（ellipse,1）
A=ellipse（k,1）;
%Amplitudechangeforthisellipsoid
asq=ellipse（k,2）^2;
%a^2
bsq=ellipse（k,3）^2;
%b^2
csq=ellipse（k,4）^2;
%c^2
x0=ellipse（k,5）;
%xoffset
y0=ellipse（k,6）;
%yoffset
z0=ellipse（k,7）;
%zoffset
phi=ellipse（k,8）*pi/180;
%firstEulerangleinradians
theta=ellipse（k,9）*pi/180;
%secondEulerangleinradians
psi=ellipse（k,10）*pi/180;
%thirdEulerangleinradians
cphi=cos（phi）;
sphi=sin（phi）;
ctheta=cos（theta）;
stheta=sin（theta）;
cpsi=cos（psi）;
spsi=sin（psi）;
%Eulerrotationmatrix
alpha=[cpsi*cphi-ctheta*sphi*spsicpsi*sphi+ctheta*cphi*spsispsi*stheta;
-spsi*cphi-ctheta*sphi*cpsi-spsi*sphi+ctheta*cphi*cpsicpsi*stheta;
stheta*sphi-stheta*cphictheta];
%rotatedellipsoidcoordinates
coordp=alpha*coord;
idx=find（（coordp（1,:
）-x0）.^2./asq+（coordp（2,:
）-y0）.^2./bsq+（coordp（3,:
）-z0）.^2./csq<
=1）;
p（idx）=p（idx）+A;
end
p=reshape（p,[nnn]）;
return;
functionout=flatten（in）
out=reshape（in,[1prod（size（in））]）;
function[e,n]=parse_inputs（varargin）
%eisthem-by-10arraywhichdefinesellipsoids
%nisthesizeofthephantombrainimage
n=128;
%Thedefaultsize
e=[];
defaults={'
shepp-logan'
'
modifiedshepp-logan'
yu-ye-wang'
fori=1:
nargin
ifischar（varargin{i}）%Lookforadefaultphantom
def=lower（varargin{i}）;
idx=strmatch（def,defaults）;
ifisempty（idx）
eid=sprintf（'
Images:
%s:
unknownPhantom'
mfilename）;
msg='
Unknowndefaultphantomselected.'
error（eid,'
%s'
msg）;
end
switchdefaults{idx}
case'
e=shepp_logan;
e=modified_shepp_logan;
e=yu_ye_wang;
elseifnumel（varargin{i}）==1
n=varargin{i};
%ascalaristheimagesize
elseifndims（varargin{i}）==2&
&
size（varargin{i},2）==10
e=varargin{i};
%userspecifiedphantom
else
eid=sprintf（'
invalidInputArgs'
msg='
Invalidinputarguments.'
error（eid,'
end
%ellipseisnotyetdefined
ifisempty（e）
e=modified_shepp_logan;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Defaultheadphantoms:
%
functione=shepp_logan
e=modified_shepp_logan;
e（:
1）=[1-.98-.02-.02.01.01.01.01.01.01];
functione=modified_shepp_logan
%ThisheadphantomisthesameastheShepp-Loganexcept
%theintensitiesarechangedtoyieldhighercontrastin
%theimage.TakenfromToft,199-200.
%
%Aabcx0y0z0phithetapsi
%-----------------------------------------------------------------
e=[1.6900.920.810000000
-.8.6624.874.7800-.01840000
-.2.1100.310.220.2200-18010
-.2.1600.410.280-.220018010
.1.2100.250.4100.35-.15000
.1.0460.046.0500.1.25000
.1.0460.046.0500-.1.25000
.1.0460.023.050-.08-.6050000
.1.0230.023.0200-.6060000
.1.0230.046.020
升级会员 