210918-放疗gamma通过率matlab代码

最新推荐文章于 2021-12-23 13:59:23 发布

xiekai1116

最新推荐文章于 2021-12-23 13:59:23 发布

阅读量1.8k

点赞数

分类专栏： dicom 文章标签： matlab

本文链接：https://blog.youkuaiyun.com/csdnxiekai/article/details/120362826

版权

dicom 专栏收录该内容

15 篇文章

订阅专栏

本文介绍了一种用于放射治疗的剂量分布分析方法——绝对剂量伽马分析，并提供了详细的MATLAB实现代码。该方法用于比较参考剂量分布图与计算剂量分布之间的差异，以评估放疗计划的准确性。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

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

gamma2.m

function [Ir,r] = gamma2(img,I,ra,dta,ps)
%绝对剂量伽马分析；img为参考剂量分布图，I为计算剂量分布，Ir为gamma因子分布，r为通过率,ps为pixelspacing
maximg=double(max(img(:))); %小数除以整数会得到错误结果
[sx,sy]=size(img);
img=double(img);I=double(I);ra=double(ra);dta=double(dta);ps=double(ps);
%Ir=zeros(512,512);
Ir=zeros(sx,sy);
detaD=ra*0.01;   %剂量差异标准3%
detad=dta;      %DTA标准3mm 
area=img>maximg*0.05;  %统计有剂量的区域，排除体表外区域  5%
%img=1/max(img(:))*img;  %相对分布
%I=1/max(I(:))*I;
%for i=1:512
    %for j=1:512
for i=1:sx
   for j=1:sy
        
        if area(i,j)==1
                      %对10个像素内的计算,ms,md行的起止处，ns,nd列的起止处
            if i-10<1
                ms=1;  
            else
                ms=i-10;
            end
            if j-10<1
                ns=1;
            else
                ns=j-10;
            end
           % if i+10>512
               % md=512;
            if i+10>sx
                md=sx;
            else
                md=i+10;
            end
            %if j+10>512
              %  nd=512;
            if j+10>sy
                nd=sy;
            else
                nd=j+10;
            end

            d=zeros(md-ms+1,nd-ns+1);
            D=zeros(md-ms+1,nd-ns+1);
            ga=zeros(md-ms+1,nd-ns+1);
            
            for k=ms:md
                for p=ns:nd
                    kk=k-ms+1;pp=p-ns+1;
                    %d(kk,pp)=sqrt((k-i)^2+(p-j)^2)*0.7828;    %pixelspacing:0.7828mm
                    d(kk,pp)=sqrt((k-i)^2+(p-j)^2)*ps;
                    %D(kk,pp)=(I(k,p)-img(i,j))/img(i,j);%global or local
                    D(kk,pp)=(I(k,p)-img(i,j))/maximg;
                    ga(kk,pp)=sqrt(d(kk,pp)^2/detad^2+D(kk,pp)^2/detaD^2);
                end
            end
          
          Ir(i,j)=min(ga(:));
        end
        
    end
end
r=length(find(Ir>1))/sum(area(:));
r=1-r;