程序流圖:

MATLAB源碼：

ClC；

clearall;

closeall;

%loadmri載入mri數(shù)據(jù)，是matlab自帶庫壓

%ph=squeeze(D);縮載入的數(shù)據(jù)D,并賦值給ph

ph=pha

ntom3d(128);

prompt={'Enterthe提示信息“輸入1到27的片的數(shù)字”

name='lnputnu彈出框名稱

哨鼬a默認(rèn)數(shù)字

R郁郁司t'&putdlg(prompt,name,1,defaulta彈出框，并得到用戶的輸入信息

nswer)%

P=squeeze(ph(:,:,str2num(cell2mat(numln將用戶的輸入信息轉(zhuǎn)換成數(shù)字，并

put))));%從ph中得到相應(yīng)的片信息P

imshow(P)展示圖片P

D=250;將D賦值為250,是從扇束頂點到旋轉(zhuǎn)

中心的像素距離.

dsensorl=2;%正實數(shù)指定扇束傳感器的間距2

F1=fanbeam(P,D,'FanSensorSpacing'.dsensorl);%通過P,D等計算扇束的數(shù)據(jù)值

dsensor2=1;%正實數(shù)指定扇束傳感器的間距1

F2=fanbeam(P,D,'FanSensorSpacing'.dsensor2);%通過P,D等計算扇束的數(shù)據(jù)值

dsensor3=%正實數(shù)指定扇束傳感器的間距0.25

[F3^sensor_pos3,fan_rotangles3]=fanbeam(P,D,...

'FanSensorSpacing',dsensor3);%通過P,D等計算扇束的數(shù)據(jù)值，并

得到扇束傳感器的位置sensor_pos3和旋轉(zhuǎn)角度fan_rot_angles3

figure,%創(chuàng)建窗口

imagesc(fan_rot_angles3,sensor_pos3,F3)%根據(jù)計算出的位置和角度展示=3的

圖片colormap(hot);

colorbar;設(shè)置色圖為hot

顯示色欄

xlabelfFanRotationAngle(degrees)')ylabel('FanSensor定義x坐標(biāo)軸

Position(degrees)')output_size=max(size(P));output_size定義y坐標(biāo)軸

Ifanl=ifanbeam(F1,D,...得到P維數(shù)的最大值，并賦值給

'FanSensorSpacing',dse

nsorl/OutputSize'.outpuLsize);

%

根據(jù)扇束投影數(shù)據(jù)F1及D等數(shù)據(jù)重

建圖像

figure,imshow(lfan1)titlef圖一');創(chuàng)建窗口，并展示圖片Ifanl

disp('圖一和原圖的性噪比為：’)；

resull=psnr1(lfari1,P);

Ifan2=ifanbeam(F2,D,...

'FanSensorSpacing',dse

nsor2,'OutputSize',output_size);

根據(jù)扇束投影數(shù)據(jù)F2及D等數(shù)據(jù)重

%

建圖像

figure,imshow(lfan2)創(chuàng)建窗口，并展示圖片Ifan2

disp('圖二和原圖的性噪比為：’)；

result=psnr1(lfan2,P);

title('圖二')；

Ifan3=ifanbeam(F3,D,...

'FanSensorSpacing',dse

nsorS.'OutputSize'.outpuLsize);根據(jù)扇束投影數(shù)據(jù)F3及D等數(shù)據(jù)重

建圖像％

figure,imshow(lfa創(chuàng)建窗口，并展示圖片Ifan3

n3)titleC圖三')；

disp（'圖三和原圖的性噪比為：'）；

result=psnr1(lfan3,P);

function[p,ellipse]=phantorr3d(varargin)%PHANT0M3DTh^ee-dimensionalan

alogueofMATLABShepp-Loganphantom%P=PHANTOM3D(DEF,N)generates

a3Dheadphantomthatcan

%beusedtotest3-Dreconstructionalgorithms.

%

%DEFisastringthatspecifiesthetypeofheadphantomtogenerate.

%Validvaluesare:%

%'Shepp-Logan'Atestimageusedwidelybyresearchersin

%tomography

%'ModifiedShepp-Logan'(default)AvariantoftheShepp-Loganpha

ntom

%inwhichthecontrastisimprovedforbetter

%visualperception.

%%NisascalarthatspecifiesthegridsizeofP.

%Ifyouomittheargument,Ndefaultsto64.

%

%P=PHANTOM3D(E,N)generatesauser-definedphantom,whereeachrow

%ofthematrixEspecifiesanellipsoidintheimage.Ehaslencolumns,

%witheachcolumncontainingadifferentparameterfortheellipsoids:

%

%Column1:Atheadditiveintensityvalueoltheellipsoid

%Column2:athelengthofthexsemi-axisoftheellipsoid

%Column3:bthelengthoftheysemi-axisoftheellipsoid

%Column4:cthelengthofthezsemi-axisoftheellipsoid

%Column5:xOthex-coordinateofthecenteroftheellipsoid

%Column6:yothey-coordinateofthecenteroftheellipsoid

%Column7:zOthez-coordinateofthecenteroftheellipsoid

%Column8:phiphiEulerangle(indegrees)(rotationaboutz-axis)

thetaEulerangle(indegrees)(rotationaboutx-

%Column9:thetaAA\

%ColumnpsipsiEulerangle(indegrees)(rotationaboutz-axis)

%

%Forpurposesofgeneratingthephantom,thedomainsforthex-,y-,and

%z-axesspan[-1,1].Columns2through7mustbespecifiedinterms%ofthisran

ge.

%

%[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

%UtahCenterforAdvaneedImagingResearch

%729ArapeenDrive

%SaltLakeCity,UT84108-1218

%

%ThiscodeisreleasedundertheGnuPublicLicense(GPL).Formoreinformation,

%see:http://www.g/copyleft/gpl.html

%

%Portionsofthiscodearebasedonphantom.m,copyrightedbytheMathworks

%

[ellipse,n]=parsejnputs(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(7);flatten(z)];

p=flatten(p);

fork=1:size(ellipse,1)

A=ellipse(k,1);%Amplitudechangeforthsellipsoid

asq=ellipse(k,2)A2;%aA2

bsq=ellipse(k,3)A2;%bA2

esq=ellipse(k,4)A2;%cA2

xOellipse(k,5);%xoffset

yoellipse(k,6);%yoffset

zOellipse(k,7);%zoffset

phi=ellipse(k,8)*pi/180;%firstEulerangleinradians

theta=ellipse(k,9)*pi/180;%secondEulerangleinradia

psi=ellipse(k,10)*pi/180;%thirdEulerangleinradia

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*cphi

ctheta];

%rotatedellipsoidcoordinates

coordp=alpha'coord;

idx=find((coordp(1,:)-x0)42./asq+(coordp(2,:)-y0)42./bsq+(coordp(3,:)-zO).A2./csq<=1);

p(idx)=p(idx)+A;

endp=reshape(p,[nnn]);

return;

functionout=flatten(in)

out=reshape(in,[1prod(size(in))]);

return;

function[e,n]=parse_inputs(varargin)

%eisthem-by-10arraywhichdefinesellipsoids

%nisthesizeofthephantombrainimage

n=128;%Thedefaultsize

e=[]；

defaults={'shepp-logan','modifiedshepp-logan','yu-ye-warg'};

fori=1:nargin

ifischar(varargin{i})%Lookforadefaultphantom

def=lower(varargin(i));

idx=strmatch(def,defaults);

ifisempty(idx)

eid=sprintf('lmages:%s:unknownPhantom',mfilename);msg='Unknowndefault

phantomselected.';error(eid,'%s',msg);

end

switchdefaults{idx}

case'shepp-logan'

e=shepplogan;

case'modifiedshepp-logan'

e=modifiedshepplogan;

case'yu-ye-wang'

e=yu_ye_wang;

end

elseifnumel(varargin{i})==1

n=varargin{i};%ascalaristheimagesize

elseifndims(varargin{i})==2&&size(varargin{i},2)==10

e=varargin{i};%userspecifiedphantom

else

eid=sprintf('lmages:%s:invalidlnputArgs'.mfilename);

msg='Invalidinputarguments.';error(eid,'%s',msg);

end

end

%ellipseisnotyetdefined

ifisempty(e)

e=modified_shepp_logan;

end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Defaultheadphantoms:%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

functione=sheppjogane=modified_shepp_logan;

c(:,1)=[1-.98.02-.02.01.01.01.01.01.01];

return;

functione=modifiedshepplogan

%

%ThisheadphantomisthesameastheShepp-Loganexcept

%theintensitiesarechangedtoyieldhighercontrastin

%theimage.TakenfromToft,199-200.

%

%AabcxOyOzOphithetapsi

%............................—

e=(1.6900.9201.810000000

-.8.6624.874.780o-.01840000

-.2.1100.310.220.2200-18010

-.2.1600.410.280-.22o018010

.1.2100.250.4100.35-.15000

.1.0460.046.0500.1.25000

.1.0460.046.050o--1.25000

.1.0460.023.050-..08-6050000

.1.0230.023.0200-.6060000

.1.0230.046.020.06--6050000];

return;

functione=yu_ye.wang

%

%YuH,YeY,WangG,Katsevich-TypeAlg