function func_conebeam_projection_Joseph(P, SOD, SDD,Nx, Ny, Nz, Pixel_size, detector_Nh,detector_Nv, detector_pixel_size, angles, filename)
% Fan beam projection function
% ------------------------------
% 算法思路:
% 1. 假设被投影物体不动,光源和探测器围绕物体逆时针旋转;
% 2. 探测器等效到原点位置,即使用虚拟探测器以便于计算;
% 3. 选定探测器通道 (hh,vv),利用射线源和该通道的坐标既可确定该射线方程(两点式形式);
% 4. 将3D物体视为由多个竖直平面组成的,具体是与x轴平行的平面还是与y轴平行的平面根据射线角度决定;
% 5. 计算射线方程与各个竖直平面的交点;
% 6. 确定交点周围四个点或八个点的图像像素标号;
% 7. 通过线性差值确定该点处的投影贡献;
% 8. 进行射线长度校正(利用角度余弦)。
% 术语:
% in-plane : 即x-o-y平面
% out-of-plane: 垂直于x-o-y平面的射线所在平面
% ------------------------------
% Inpit paramaters:
% P : 3D投影图像(即S_L头模型产生的三维图像矩阵),坐标顺序:P(nx,ny,nz)
% SOD : source object distance
% SDD : source detector distance
% Nx,Ny ,Nz:三维图像的大小
% Pixel_size :一个像素大小(即体素的一个面)
% detector_Nh : 探测器的水平通道个数
% detector_Nv : 探测器的竖直通道个数
% detector_pixel_size
% angles : 投影角度矢量 in degree
% ------------------------------
% Output paramaters:
% g : projection data (Nh*Nv), output by file
% g = zeros(detector_Nh, detector_Nv); % 探测器大小为 Nh(宽,即列数) * Nv(高,即行数)
% 由于matlab是列优先排布,因此g的实际排布形式如下&