模拟的r射线，不用考虑散射了。
function []=spheresource()h=input('source to collimator:');H=input('source to detector:');r=input('hole diameter:');d1=input('pinhole1 parameter(x y axis):');d2=input('pinhole2 parameter:');d3=input('pinhole3 parameter:');
d4=input('pinhole4 parameter:');%h is distanse source to collimator,H is source to detector%dx is center of the hole to y axis%dy is center of the hole to x axis%spheresource(5,20,3,3,4)%spheresource(5,20,-3,3,4)%spheresource(5,20,2.5 ,2.5,-2.5 ,2.5,2.5,-2.5,-2.5,-2.5,4)
N=80;%6可调分辨率I64=zeros(N,N); %预定义平面164的灰度值为O[m,n]=meshgrid(linspace(-N/2,N/2-1,N));%确定坐标系及坐标原点的位置
% h=1;H=2;% dx=5;% r=5;
r1=r*H/h; %圆半径大小控制变量(单位pixel)a1=d1(1)*H/h; %圆心位置控制变量b1=d1(2)*H/h; %圆心位置控制变量D1=((m+a1).^2+(n+b1).^2).^(1/2); %圆函数关系式i=find(D1<=r1); %单下标寻址实现语句，find()函数足条件D<-r的像素点的单下标值I64(i)=0.2; %按要求为像素点赋值
r2=r*H/h; %圆半径大小控制变量(单位pixel)a2=d2(1)*H/h; %圆心位置控制变量b2=d2(2)*H/h; %圆心位置控制变量D2=((m+a2).^2+(n+b2).^2).^(1/2); %圆函数关系式 i=find(D2<=r2);%单下标寻址实现语句，find()函数足条件D<-r的像素点的单下标值
I64(i)=0.2+I64(i); %按要求为像素点赋值
r3=r*H/h; %圆半径大小控制变量(单位pixel)a3=d3(1)*H/h; %圆心位置控制变量b3=d3(2)*H/h; %圆心位置控制变量D3=((m+a3).^2+(n+b3).^2).^(1/2); %圆函数关系式 i=find(D3<=r3);%单下标寻址实现语句，find()函数足条件D<-r的像素点的单下标值
I64(i)=0.2+I64(i); %按要求为像素点赋值
r4=r*H/h; %圆半径大小控制变量(单位pixel)a4=d4(1)*H/h; %圆心位置控制变量b4=d4(2)*H/h; %圆心位置控制变量D4=((m+a4).^2+(n+b4).^2).^(1/2); %圆函数关系式 i=find(D4<=r4);%单下标寻址实现语句，find()函数足条件D<-r的像素点的单下标值
I64(i)=0.2+I64(i); %按要求为像素点赋值
J=mat2gray(I64,[0,1]);%matrix to gray scale imageimshow(J);len=length(find(I64(i)==0.8))%find the gray scale equals to 0.8(4order overlap)overlap=len/(N*N)%percentage
>> spheresource()source to collimator:3source to detector:10hole diameter:1pinhole1 parameter(x y axis):[2 2]pinhole2 parameter:[-2 2]pinhole3 parameter:[2 -2]pinhole4 parameter:[-2 -2]
len =
0
overlap =
0
