210918-放疗gamma通过率matlab代码

Good

引文：放疗中常用通过率，一种是使用matrixx软件(比较平面)，ArcCheck软件(比较探头所在平面)，另一种是使用代码，这里仍然来自Niubi高，做此记录，不然下次就找不到了。

gamma2.m

1. function [Ir,r] = gamma2(img,I,ra,dta,ps)
2. %绝对剂量伽马分析；img为参考剂量分布图，I为计算剂量分布，Ir为gamma因子分布，r为通过率,ps为pixelspacing
3. maximg=double(max(img(:))); %小数除以整数会得到错误结果
4. [sx,sy]=size(img);
5. img=double(img);I=double(I);ra=double(ra);dta=double(dta);ps=double(ps);
6. %Ir=zeros(512,512);
7. Ir=zeros(sx,sy);
8. detaD=ra*0.01;   %剂量差异标准3%
9. detad=dta;      %DTA标准3mm
10. area=img>maximg*0.05;  %统计有剂量的区域，排除体表外区域  5%
11. %img=1/max(img(:))*img;  %相对分布
12. %I=1/max(I(:))*I;
13. %for i=1:512
14.     %for j=1:512
15. for i=1:sx
16.    for j=1:sy
17.         
18.         if area(i,j)==1
19.                       %对10个像素内的计算,ms,md行的起止处，ns,nd列的起止处
20.             if i-10<1
21.                 ms=1;  
22.             else
23.                 ms=i-10;
24.             end
25.             if j-10<1
26.                 ns=1;
27.             else
28.                 ns=j-10;
29.             end
30.            % if i+10>512
31.                % md=512;
32.             if i+10>sx
33.                 md=sx;
34.             else
35.                 md=i+10;
36.             end
37.             %if j+10>512
38.               %  nd=512;
39.             if j+10>sy
40.                 nd=sy;
41.             else
42.                 nd=j+10;
43.             end
45.             d=zeros(md-ms+1,nd-ns+1);
46.             D=zeros(md-ms+1,nd-ns+1);
47.             ga=zeros(md-ms+1,nd-ns+1);
48.             
49.             for k=ms:md
50.                 for p=ns:nd
51.                     kk=k-ms+1;pp=p-ns+1;
52.                     %d(kk,pp)=sqrt((k-i)^2+(p-j)^2)*0.7828;    %pixelspacing:0.7828mm
53.                     d(kk,pp)=sqrt((k-i)^2+(p-j)^2)*ps;
54.                     %D(kk,pp)=(I(k,p)-img(i,j))/img(i,j);%global or local
55.                     D(kk,pp)=(I(k,p)-img(i,j))/maximg;
56.                     ga(kk,pp)=sqrt(d(kk,pp)^2/detad^2+D(kk,pp)^2/detaD^2);
57.                 end
58.             end
59.          
60.           Ir(i,j)=min(ga(:));
61.         end
62.         
63.     end
64. end
65. r=length(find(Ir>1))/sum(area(:));
66. r=1-r;

复制代码

来源：https://blog.csdn.net/csdnxiekai/article/details/120362826
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！