本發(fā)明屬于計(jì)算機(jī)信息技術(shù)在核技術(shù)領(lǐng)域的應(yīng)用，涉及到腫瘤放射治療技術(shù)，具體是基于gpu蒙特卡洛算法的磁場下質(zhì)子和重離子劑量計(jì)算方法，可用于傳統(tǒng)以及磁共振實(shí)時(shí)引導(dǎo)質(zhì)子與重離子放射治療計(jì)劃系統(tǒng)的劑量計(jì)算。

背景技術(shù)：

放射治療的目的是在最大程度殺滅腫瘤的同時(shí)保護(hù)周圍危及器官。為了實(shí)現(xiàn)這個(gè)目的，可從影像技術(shù)與選用合適的輻射源類型的兩個(gè)方面考慮。

一方面醫(yī)生需要借助影像技術(shù)清晰地定位腫瘤和周圍危及器官的位置。目前大多采用ct確定位置，但是ct對軟組織成像分辨率不高，且在ct掃描過程中，患者需要承擔(dān)放射治療以外的電離輻射。另外當(dāng)采用ct定位運(yùn)動器官位置時(shí)，一般需要在患者體內(nèi)植入標(biāo)記物，而很多患者都不愿意或者不適合接受外部植入標(biāo)記物。核磁共振成像(mri)能實(shí)時(shí)跟蹤腫瘤位置的變化，無電離輻射、不需外部植入標(biāo)記物，所以國際上提出用磁共振實(shí)時(shí)引導(dǎo)放射治療，在治療前、治療中通過mri來引導(dǎo)整個(gè)放射治療的過程。

另一方面，目前臨床上選用治療的輻射源類型多為光子和電子。電子主要適用于淺表腫瘤治療。光子治療時(shí)，對于較深處的腫瘤，在殺滅腫瘤的同時(shí)也極大地?fù)p傷了周圍健康的組織器官。質(zhì)子或重離子放射治療是目前最為先進(jìn)的放療技術(shù)之一，它充分利用了質(zhì)子或重離子劑量分布的布拉格峰特性，將絕大部分輻射劑量沉積在腫瘤靶區(qū)，同時(shí)極大地降低了腫瘤后部正常細(xì)胞或敏感器官受到的輻射損傷。

充分考慮以上兩個(gè)方面，若采用磁共振實(shí)時(shí)引導(dǎo)的質(zhì)子放射治療或者磁共振實(shí)時(shí)引導(dǎo)的重離子放射治療，不僅能給出腫瘤和周圍危及器官，包括運(yùn)動器官的清晰圖像，同時(shí)也能極大地降低對正常組織的損傷。但當(dāng)質(zhì)子和重離子處于mri磁場中時(shí)，由于其帶電，磁場將影響質(zhì)子或重離子的運(yùn)動軌跡，從而改變其在患者體內(nèi)的劑量分布，因此在設(shè)計(jì)相應(yīng)治療計(jì)劃的同時(shí)需要考慮磁場這一因素。對于某些運(yùn)動器官的腫瘤或是兒童患者，若能采用磁共振實(shí)時(shí)引導(dǎo)的質(zhì)子和重離子精確放射治療，將能進(jìn)一步提升放射治療的效果，改善病人生存質(zhì)量。

目前大多醫(yī)院在實(shí)施臨床放射治療時(shí)，治療計(jì)劃系統(tǒng)中的輻射劑量計(jì)算都是采用解析算法來得到的，主要就是筆形束卷積疊加算法。但對于一些密度差別較大的非均勻組織，筆形束算法計(jì)算得到的劑量結(jié)果并不準(zhǔn)確，這可能會導(dǎo)致靶區(qū)欠劑量而周圍危及器官劑量過高的現(xiàn)象。而對于磁場作用下的質(zhì)子和重離子輻射劑量，目前還無法根據(jù)現(xiàn)有的筆形束算法準(zhǔn)確獲得質(zhì)子和重離子劑量的分布情況。

蒙特卡洛方法一直被公認(rèn)為輻射劑量計(jì)算領(lǐng)域內(nèi)的“黃金準(zhǔn)則”，可以實(shí)現(xiàn)完全真實(shí)地模擬磁共振實(shí)時(shí)引導(dǎo)質(zhì)子和重離子放射治療的情況，計(jì)算出質(zhì)子和重離子對腫瘤靶區(qū)以及周圍危及器官的輻射劑量。但使用傳統(tǒng)的cpu串行蒙特卡洛程序計(jì)算劑量花費(fèi)時(shí)間太長，故無法真正應(yīng)用于臨床。隨著計(jì)算機(jī)硬件的發(fā)展，基于kepler構(gòu)架和maxwell構(gòu)架的gpu在當(dāng)前的科學(xué)計(jì)算中得到了廣泛的應(yīng)用。但將已有的cpu大型蒙特卡洛程序如geant4、mcnp移植到gpu上，需花費(fèi)大量的時(shí)間來改寫代碼，并且將gpu上有限的內(nèi)存分配到幾千個(gè)線程，對于蒙特卡洛方法而言，這是遠(yuǎn)遠(yuǎn)不足的。另外，基于kepler構(gòu)架和maxwell構(gòu)架的gpu缺乏固有雙精度原子加法運(yùn)算的硬件支持。傳統(tǒng)的解決方法是使用nvidia的“比較-交換”算法從軟件上效仿雙精度原子加法運(yùn)算。這個(gè)算法的特點(diǎn)是步驟不定長，gpu的線程競爭越激烈，所需步驟越多。當(dāng)前gpu硬件上采用“單指令-多線程”技術(shù)，一個(gè)線程包內(nèi)的32個(gè)線程在時(shí)間上以“鎖步”(lockstep)方式執(zhí)行指令。這樣導(dǎo)致的結(jié)果是：如果32個(gè)線程有某些需要更新同一處全局內(nèi)存，則會產(chǎn)生激烈的線程競爭，“比較-交換”算法所需步驟急劇增加，運(yùn)算時(shí)間大幅增長。因此，將gpu與蒙特卡洛程序相結(jié)合，若不解決這一系列問題，將無法快速準(zhǔn)確的計(jì)算得到劑量結(jié)果。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明是為了解決上述現(xiàn)有技術(shù)存在的不足之處，提出一種基于gpu蒙特卡洛算法的磁場下質(zhì)子和重離子劑量計(jì)算方法，以期能快速且準(zhǔn)確地計(jì)算出磁場下質(zhì)子和重離子與物質(zhì)相互作用的輻射劑量，可用于傳統(tǒng)的和磁共振實(shí)時(shí)引導(dǎo)的質(zhì)子與重離子放射治療計(jì)劃的劑量計(jì)算，進(jìn)而提高質(zhì)子和重離子劑量計(jì)算的準(zhǔn)確性與計(jì)算速度，改善放射治療的效果。

本發(fā)明為解決技術(shù)問題采用如下技術(shù)方案：

本發(fā)明一種基于gpu蒙特卡洛算法的磁場下質(zhì)子和重離子劑量計(jì)算方法的特點(diǎn)是按如下步驟進(jìn)行：

步驟1：采集數(shù)據(jù)；

步驟1.1、獲取放射治療質(zhì)子或重離子加速器的束流源信息其中，e表示源能量，表示源位置，表示源發(fā)射方向，σ是高斯分布的標(biāo)準(zhǔn)差；

步驟1.2、獲取人體解剖結(jié)構(gòu)的圖像數(shù)據(jù)并重建人體模型；獲取核磁共振儀的磁場強(qiáng)度數(shù)據(jù)

步驟1.3、獲取質(zhì)子和重離子分別與物質(zhì)發(fā)生反應(yīng)的核數(shù)據(jù)并進(jìn)行處理，得到所述核數(shù)據(jù)的宏觀截面數(shù)據(jù)∑，對所述宏觀截面數(shù)據(jù)∑按照能量的高低進(jìn)行降序排序，得到排序后的質(zhì)子宏觀截面數(shù)據(jù)∑^p、重離子宏觀截面數(shù)據(jù)∑^ion、次級粒子宏觀截面數(shù)據(jù)∑^psec＝{∑^e,∑^d,∑^α}；∑^e為次級電子宏觀截面數(shù)據(jù)；∑^d為氘核宏觀截面數(shù)據(jù)；∑^α為α粒子宏觀截面數(shù)據(jù)；

步驟2、確定gpu的最優(yōu)線程數(shù)和輸運(yùn)任務(wù)的批次；

步驟2.1、利用runtimeattribute程序接口獲得gpu中每個(gè)線程所需寄存器的數(shù)目r；則gpu中每個(gè)流多處理器工作在滿載狀態(tài)的最小線程個(gè)數(shù)為r表示每個(gè)流多處理器上的寄存器個(gè)數(shù)；從而得到gpu工作在滿載狀態(tài)所需線程總數(shù)為t＝mt，m表示gpu中流多處理器的個(gè)數(shù)；

步驟2.2、設(shè)置放射源粒子的數(shù)目為n，并將n個(gè)放射源粒子的輸運(yùn)任務(wù)劃分為t個(gè)批次，使得每個(gè)批次上以串行地方式待計(jì)算放射源粒子的個(gè)數(shù)為

步驟3、利用蒙特卡洛算法計(jì)算每個(gè)批次在磁場作用下的初級質(zhì)子和重離子輻射劑量；

步驟3.1、定義每一批抽取初級放射源粒子的次數(shù)為w，并初始化w＝1；

定義次級粒子的次數(shù)為w′，并初始化w′＝1；

步驟3.2、定義第w次抽取放射源粒子時(shí)的輸運(yùn)次數(shù)為u，并初始化u＝0；

定義第w′次抽取次級粒子時(shí)的輸運(yùn)次數(shù)為u′，并初始化u′＝0；

步驟3.3、利用隨機(jī)數(shù)生成器從放射源信息s中第w次抽取第w個(gè)放射源粒子sw；所述第w個(gè)放射源粒子sw第u次輸運(yùn)的狀態(tài)為：

步驟3.4、判斷第w次抽取第w個(gè)放射源粒子的類型，若為質(zhì)子記為則基于排序后的質(zhì)子宏觀截面數(shù)據(jù)σ^p抽取第u+1次輸運(yùn)時(shí)對應(yīng)粒子類型的運(yùn)動步長和運(yùn)動方向若為重離子記為則基于排序后的重離子宏觀截面數(shù)據(jù)σ^ion抽取第u+1次輸運(yùn)時(shí)對應(yīng)粒子類型的運(yùn)動步長和運(yùn)動方向

步驟3.5、判斷所述放射源粒子是否處于人體模型的磁場區(qū)，若是，先將第w個(gè)放射源粒子沿運(yùn)動方向移動的距離，再將第w個(gè)放射源粒子沿著式(1)修正的運(yùn)動方向移動的距離；否則，仍然采用所述運(yùn)動方向

式(1)中，為質(zhì)子或重離子第u+1次輸運(yùn)時(shí)修正后的新方向；norm{}為歸一化算符；du+1為第u+1次對第w個(gè)質(zhì)子或重離子進(jìn)行抽樣所得到的步長；q為質(zhì)子或重離子的電荷數(shù)；c為真空中光速，m為質(zhì)子或重離子質(zhì)量，為第w個(gè)質(zhì)子或重離子在第u次輸運(yùn)時(shí)的能量；為質(zhì)子或重離子第u+1次輸運(yùn)時(shí)未修正的運(yùn)動方向；為第w個(gè)質(zhì)子或重離子在人體模型中第u+1次輸運(yùn)時(shí)所處位置的磁場強(qiáng)度；

步驟3.6、基于排序后的質(zhì)子宏觀截面數(shù)據(jù)σ^p和重離子宏觀截面數(shù)據(jù)σ^ion、所述第w個(gè)放射源粒子sw第u次輸運(yùn)的狀態(tài)運(yùn)動步長、運(yùn)動方向或修正的運(yùn)動方向，對第w個(gè)放射源粒子sw的反應(yīng)類型進(jìn)行抽樣，得到第w個(gè)放射源粒子sw在人體模型中進(jìn)行第u+1次輸運(yùn)的狀態(tài)以及第u次輸運(yùn)過程產(chǎn)生的個(gè)次級粒子并存儲；

步驟3.7、計(jì)算第w次抽取的第w個(gè)放射源粒子sw在人體模型中進(jìn)行第u+1次輸運(yùn)時(shí)沉積的劑量

步驟3.8、將u+1賦值給u，并返回步驟3.4執(zhí)行，直到第w個(gè)放射源粒子sw第u次輸運(yùn)的能量低于截止能量或第w個(gè)放射源粒子sw超出人體模型的邊界為止，記第w個(gè)放射源粒子sw的輸運(yùn)總次數(shù)為u，從而統(tǒng)計(jì)得到第w次抽取的第w個(gè)放射源粒子sw在人體模型中沉積的劑量dosew；

步驟3.9、將w+1賦值給w，并返回步驟3.1執(zhí)行，直到w＞n為止；從而統(tǒng)計(jì)得到每個(gè)批次n次抽取的n個(gè)放射源粒子在人體模型中沉積的劑量dose；

步驟3.10、統(tǒng)計(jì)得到所有t個(gè)批次的n個(gè)放射源粒子輸運(yùn)u次產(chǎn)生的次級粒子的數(shù)目為并將n′個(gè)次級粒子的輸運(yùn)任務(wù)劃分為t個(gè)批次，使得每個(gè)批次上以串行地方式待計(jì)算次級粒子的個(gè)數(shù)為

步驟4、利用蒙特卡洛算法計(jì)算每個(gè)批次在磁場作用下的次級粒子的輻射劑量；

步驟4.1、判斷n′個(gè)次級粒子是否是在質(zhì)子輸運(yùn)過程產(chǎn)生的，若是，則表示n′個(gè)次級粒子是由次級質(zhì)子、次級電子、氘核和α粒子組成，并執(zhí)行步驟4.2；否則，表示n′個(gè)次級粒子是由次級¹²c粒子、α粒子、次級質(zhì)子和次級¹¹c粒子組成，并執(zhí)行步驟4.3；

步驟4.2、對次級質(zhì)子、次級電子、氘核和α粒子分別按如下方式進(jìn)行輸運(yùn)：

若為次級質(zhì)子，則按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為次級電子，則基于次級電子宏觀截面數(shù)據(jù)∑^e，按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為氘核，則基于氘核宏觀截面數(shù)據(jù)∑^d，按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為α粒子，則基于α粒子宏觀截面數(shù)據(jù)∑^α，按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

步驟4.3、對次級¹²c粒子、α粒子、次級質(zhì)子和次級¹¹c粒子分別按如下方式進(jìn)行輸運(yùn)：

若為次級¹²c粒子，則按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為α粒子，則基于α粒子宏觀截面數(shù)據(jù)∑^α，按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為次級質(zhì)子，則按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為次級¹¹c粒子，則按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

步驟5、基于gpu快速原子加法統(tǒng)計(jì)劑量結(jié)果：

步驟5.1、將t個(gè)批次所獲得的放射源粒子和次級粒子的沉積劑量以三維矩陣的形式存入gpu的全局內(nèi)存中，在所述三維矩陣中的任意元素記為dosei；dosei表示在三維空間中第i個(gè)位置上的劑量；

步驟5.2、獲取gpu中同一個(gè)線程包內(nèi)的有效線程，并將所要更新的全局內(nèi)存地址相同的有效線程從所述同一線程包中篩選出來；

步驟5.3、將篩選出來的有效線程所對應(yīng)的全局內(nèi)存中所存儲的劑量進(jìn)行累加，得到的劑量結(jié)果存入序號最小的有效線程中；

步驟5.4、利用“卡漢求和”以及“比較-交換”算法將所有屬于同一線程包內(nèi)序號最小的有效線程中的劑量結(jié)果累加到gpu的全局內(nèi)存中，從而得到人體模型中放射源粒子和次級粒子的總劑量；

步驟5.5、將所述總劑量除以放射源粒子的總數(shù)n，從而得到歸一化的劑量結(jié)果，并將其返回到cpu內(nèi)存。

與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果在于：

1、本發(fā)明利用圖像數(shù)據(jù)構(gòu)建人體模型，并調(diào)用cudaruntimeapi里面的函數(shù)cudafuncattributes()獲得gpu中每個(gè)線程所需寄存器的數(shù)目r，進(jìn)而計(jì)算得到gpu的最優(yōu)線程分配方案；將n個(gè)放射源粒子的模擬任務(wù)在gpu上劃分為t個(gè)批次，每個(gè)批次同步執(zhí)行；基于排序后的宏觀截面數(shù)據(jù)σ，利用蒙特卡洛算法計(jì)算每個(gè)批次在磁場作用下的質(zhì)子和重離子輻射劑量，對磁場作用下的帶電粒子，經(jīng)修正運(yùn)動方向后再輸運(yùn)；最后采用一種新的gpu快速原子加法統(tǒng)計(jì)劑量結(jié)果，從而解決了目前該領(lǐng)域內(nèi)4個(gè)主要難題，包括：1.現(xiàn)有治療計(jì)劃系統(tǒng)中采用筆形束卷積疊加的解析算法無法準(zhǔn)確計(jì)算出質(zhì)子和重離子在人體模型非均勻介質(zhì)中沉積劑量的問題；2.現(xiàn)有治療計(jì)劃系統(tǒng)解析算法無法考慮磁共振實(shí)時(shí)引導(dǎo)放射治療時(shí)磁場對劑量的影響，從而導(dǎo)致計(jì)算結(jié)果不準(zhǔn)確；3.現(xiàn)有的傳統(tǒng)蒙特卡洛程序，基于cpu硬件平臺，計(jì)算速度非常慢；4.基于kepler構(gòu)架和maxwell構(gòu)架的gpu在處理大量數(shù)據(jù)求和時(shí)缺乏固有雙精度原子加法運(yùn)算硬件支持的問題。本發(fā)明解決了上述的一系列問題，可用于治療計(jì)劃系統(tǒng)中質(zhì)子和重離子輻射劑量的計(jì)算，包括在磁場作用下質(zhì)子和重離子的蒙特卡洛輸運(yùn)模擬，可改進(jìn)目前臨床治療計(jì)劃系統(tǒng)中的解析法的固有缺陷，提高了劑量計(jì)算的準(zhǔn)確性，改善質(zhì)子和重離子放射治療的效果，進(jìn)而提高了患者術(shù)后生存質(zhì)量。

2、本發(fā)明中利用gpu的最優(yōu)線程分配法完成n個(gè)放射源粒子的模擬任務(wù)，該方法利用cudaruntimeapi里面的函數(shù)cudafuncattributes()獲得gpu中每個(gè)線程所需寄存器的數(shù)目r，通過r進(jìn)一步獲得gpu中每個(gè)流多處理器工作在滿載狀態(tài)的最小線程個(gè)數(shù)t，最后計(jì)算出gpu工作在滿載狀態(tài)時(shí)所需線程總數(shù)t，解決了gpu上線程不能充分合理利用的問題，從而能夠最大化gpu的并行能力，提高了蒙特卡洛輸運(yùn)模擬的計(jì)算速度。

3、本發(fā)明在gpu上輸運(yùn)粒子之前，對宏觀截面數(shù)據(jù)σ按照能量排序，從而提高了實(shí)際模擬粒子輸運(yùn)任務(wù)時(shí)查找宏觀截面數(shù)據(jù)的速度，進(jìn)而提高了模擬的速度。

4、本發(fā)明在gpu上模擬粒子的輸運(yùn)任務(wù)時(shí)，利用蒙特卡洛算法計(jì)算每個(gè)批次在磁場作用下的質(zhì)子和重離子輻射劑量，尤其針對磁場會導(dǎo)致劑量的變化，在輸運(yùn)過程中做了初級帶電粒子和次級帶電粒子運(yùn)動方向的修正，優(yōu)化了原始的蒙特卡洛方法，使得新算法適用范圍更廣；從而解決了傳統(tǒng)蒙特卡洛方法無法處理磁場存在時(shí)劑量計(jì)算的問題。將本發(fā)明應(yīng)用于磁共振實(shí)時(shí)引導(dǎo)質(zhì)子和重離子的治療計(jì)劃系統(tǒng)中的劑量計(jì)算，大幅提高了質(zhì)子和重離子輻射劑量計(jì)算的準(zhǔn)確性與計(jì)算速度。

5、本方法中提出一種新的gpu快速原子加法統(tǒng)計(jì)劑量結(jié)果，采用卡漢求和算法、nvidia的韋斯福算法和nvidia的“比較-交換”算法從軟件上效仿雙精度原子加法運(yùn)算，解決了基于kepler構(gòu)架和maxwell構(gòu)架的gpu缺乏固有雙精度原子加法運(yùn)算的硬件支持問題，并且可快速準(zhǔn)確實(shí)現(xiàn)任意類型數(shù)據(jù)的計(jì)算，還減少了線程包內(nèi)部的線程競爭，提升了單精度原子加法運(yùn)算的精度和雙精度原子加法運(yùn)算的速度，因而提高了統(tǒng)計(jì)劑量結(jié)果的準(zhǔn)確性與速度。

附圖說明

圖1為本發(fā)明的整體流程圖。

具體實(shí)施方式

本實(shí)施例中，一種基于gpu蒙特卡洛算法的磁場下質(zhì)子和重離子劑量計(jì)算方法，是應(yīng)用于質(zhì)子和重離子放射治療領(lǐng)域，即質(zhì)子束與重離子束通過與人體模型的相互作用，將能量沉積在人體模型特定位置的實(shí)驗(yàn)環(huán)境中，具體的說，如圖1所示，按如下步驟進(jìn)行：

步驟1：采集數(shù)據(jù)；

步驟1.1、獲取放射治療質(zhì)子或重離子加速器的束流源信息其中，e表示源能量，表示源位置，表示源發(fā)射方向，σ是高斯分布的標(biāo)準(zhǔn)差；

步驟1.2、獲取ct機(jī)或核磁共振儀mri或正電子發(fā)射型計(jì)算機(jī)斷層顯像機(jī)器pet的人體解剖結(jié)構(gòu)的圖像數(shù)據(jù)并重建人體模型；獲取核磁共振儀的磁場強(qiáng)度數(shù)據(jù)本實(shí)施例中采集的圖像數(shù)據(jù)支持更多種格式，解決了現(xiàn)有技術(shù)中無法通過ct機(jī)獲取人體模型中軟組織、運(yùn)動器官準(zhǔn)確位置的問題，并解決了經(jīng)ct掃描獲得人體模型數(shù)據(jù)時(shí)會受到放射源帶來的輻射劑量問題，支持mri和pet格式的數(shù)據(jù)提高了在人體模型中所包含的腫瘤區(qū)域的界限精準(zhǔn)度，減少了癌癥患者為獲得人體模型數(shù)據(jù)需要進(jìn)行體內(nèi)植入標(biāo)記物的痛苦；

步驟1.3、獲取質(zhì)子和重離子分別與物質(zhì)發(fā)生反應(yīng)的核數(shù)據(jù)并進(jìn)行處理，得到所述核數(shù)據(jù)的宏觀截面數(shù)據(jù)∑，對所述宏觀截面數(shù)據(jù)∑按照能量的高低進(jìn)行降序排序，得到排序后的質(zhì)子宏觀截面數(shù)據(jù)∑^p、重離子宏觀截面數(shù)據(jù)∑^ion、次級粒子宏觀截面數(shù)據(jù)∑^psec＝{∑^e,∑^d,∑^α}；∑^e為次級電子宏觀截面數(shù)據(jù)；∑^d為氘核宏觀截面數(shù)據(jù)；∑^α為α粒子宏觀截面數(shù)據(jù)；宏觀截面數(shù)據(jù)按照能量順序存儲，便于后面反應(yīng)類型的抽取。

步驟2、確定gpu的最優(yōu)線程數(shù)和輸運(yùn)任務(wù)的批次；

步驟2.1、利用runtimeattribute程序接口獲得gpu中每個(gè)線程所需寄存器的數(shù)目r；則gpu中每個(gè)流多處理器工作在滿載狀態(tài)的最小線程個(gè)數(shù)為r表示每個(gè)流多處理器上的寄存器個(gè)數(shù)；從而得到gpu工作在滿載狀態(tài)所需線程總數(shù)為t＝mt，m表示gpu中流多處理器的個(gè)數(shù)；

步驟2.2、設(shè)置放射源粒子的數(shù)目為n，并將n個(gè)放射源粒子的輸運(yùn)任務(wù)劃分為t個(gè)批次，則在每個(gè)批次上以串行地方式待計(jì)算放射源粒子的個(gè)數(shù)為

步驟3、利用蒙特卡洛算法計(jì)算每個(gè)批次在磁場作用下的初級質(zhì)子和重離子輻射劑量；

步驟3.1、定義每一批抽取初級放射源粒子的次數(shù)為w，并初始化w＝1；

定義次級粒子的次數(shù)為w′，并初始化w′＝1；

步驟3.2、定義第w次抽取放射源粒子時(shí)的輸運(yùn)次數(shù)為u，并初始化u＝00；

定義第w′次抽取次級粒子時(shí)的輸運(yùn)次數(shù)為u′，并初始化u′＝0；

步驟3.3、利用隨機(jī)數(shù)生成器從放射源信息s中第w次抽取第w個(gè)放射源粒子sw；所述第w個(gè)放射源粒子sw第u次輸運(yùn)的狀態(tài)為：

步驟3.4、判斷第w次抽取第w個(gè)放射源粒子的類型，若為質(zhì)子記為則基于排序后的質(zhì)子宏觀截面數(shù)據(jù)∑^p抽取第u+1次輸運(yùn)時(shí)對應(yīng)粒子類型的運(yùn)動步長和運(yùn)動方向若為重離子記為則基于排序后的重離子宏觀截面數(shù)據(jù)∑^ion抽取第u+1次輸運(yùn)時(shí)對應(yīng)粒子類型的運(yùn)動步長和運(yùn)動方向

步驟3.5、判斷所述放射源粒子是否處于人體模型的磁場區(qū)，若是，先將第w個(gè)放射源粒子沿運(yùn)動方向移動的距離，再將第w個(gè)放射源粒子沿著式(1)修正的運(yùn)動方向移動的距離；否則，仍然采用所述運(yùn)動方向

式(1)中，為質(zhì)子或重離子第u+1次輸運(yùn)時(shí)修正后的新方向；norm{}為歸一化算符；du+1為第u+1次對第w個(gè)質(zhì)子或重離子進(jìn)行抽樣所得到的步長；q為質(zhì)子或重離子的電荷數(shù)；c為真空中光速，m為質(zhì)子或重離子質(zhì)量，為第w個(gè)質(zhì)子或重離子在第u次輸運(yùn)時(shí)的能量；為質(zhì)子或重離子第u+1次輸運(yùn)時(shí)未修正的運(yùn)動方向；為第w個(gè)質(zhì)子或重離子在人體模型中第u+1次輸運(yùn)時(shí)所處位置的磁場強(qiáng)度；

根據(jù)原先的蒙特卡洛方法并不能計(jì)算磁場下的質(zhì)子和重離子劑量問題，本發(fā)明考慮了磁場對質(zhì)子、重離子的影響，修正了質(zhì)子和重離子的運(yùn)動方向，進(jìn)而獲得較為準(zhǔn)確的磁場下質(zhì)子、重離子的運(yùn)動軌跡，有助于優(yōu)化原先的蒙特卡洛方法，拓寬它的應(yīng)用范圍。本算法可用于核磁實(shí)時(shí)引導(dǎo)質(zhì)子放射治療或核磁實(shí)時(shí)引導(dǎo)重離子放射治療的計(jì)劃系統(tǒng)的劑量計(jì)算。

步驟3.6、基于排序后的質(zhì)子宏觀截面數(shù)據(jù)∑^p和重離子宏觀截面數(shù)據(jù)∑^ion、所述第w個(gè)放射源粒子sw第u次輸運(yùn)的狀態(tài)運(yùn)動步長、運(yùn)動方向或修正的運(yùn)動方向，對第w個(gè)放射源粒子sw的反應(yīng)類型進(jìn)行抽樣，得到第w個(gè)放射源粒子sw在人體模型中進(jìn)行第u+1次輸運(yùn)的狀態(tài)以及第u次輸運(yùn)過程產(chǎn)生的個(gè)次級粒子并存儲；

步驟3.7、計(jì)算第w次抽取的第w個(gè)放射源粒子sw在人體模型中進(jìn)行第u+1次輸運(yùn)時(shí)沉積的劑量

步驟3.8、將u+1賦值給u，并返回步驟3.4執(zhí)行，直到第w個(gè)放射源粒子sw第u次輸運(yùn)的能量低于截止能量或第w個(gè)放射源粒子sw超出人體模型的邊界為止，記第w個(gè)放射源粒子sw的輸運(yùn)總次數(shù)為u；從而統(tǒng)計(jì)得到第w次抽取的第w個(gè)放射源粒子sw在人體模型中沉積的劑量dosew；

步驟3.9、將w+1賦值給w，并返回步驟3.1執(zhí)行，直到w＞n為止；從而統(tǒng)計(jì)得到每個(gè)批次n次抽取的n個(gè)放射源粒子在人體模型中沉積的劑量dose；

步驟3.10、統(tǒng)計(jì)得到所有t個(gè)批次的n個(gè)放射源粒子輸運(yùn)u次產(chǎn)生的次級粒子的數(shù)目為并將n′個(gè)次級粒子的輸運(yùn)任務(wù)劃分為t個(gè)批次，使得每個(gè)批次上以串行地方式待計(jì)算次級粒子的個(gè)數(shù)為

步驟4、利用蒙特卡洛算法計(jì)算每個(gè)批次在磁場作用下的次級粒子的輻射劑量；

步驟4.1、判斷n′個(gè)次級粒子是否是在質(zhì)子輸運(yùn)過程產(chǎn)生的，若是，則表示n′個(gè)次級粒子是由次級質(zhì)子、次級電子、氘核和α粒子組成，并執(zhí)行步驟4.2；否則，表示n′個(gè)次級粒子是由次級¹²c粒子、α粒子、次級質(zhì)子和次級¹¹c粒子組成，并執(zhí)行步驟4.3；

步驟4.2、對次級質(zhì)子、次級電子、氘核和α粒子分別按如下方式進(jìn)行輸運(yùn)：

若為次級質(zhì)子，則按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為次級電子，則基于次級電子宏觀截面數(shù)據(jù)∑^e，按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為氘核，則基于氘核宏觀截面數(shù)據(jù)∑^d，按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為α粒子，則基于α粒子宏觀截面數(shù)據(jù)∑^α，按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

步驟4.3、對次級¹²c粒子、α粒子、次級質(zhì)子和次級¹¹c粒子分別按如下方式進(jìn)行輸運(yùn)：

若為次級¹²c粒子，則按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為α粒子，則基于α粒子宏觀截面數(shù)據(jù)∑^α，按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為次級質(zhì)子，則按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

若為次級¹¹c粒子，則按照步驟3.4-步驟3.9進(jìn)行輸運(yùn)；

質(zhì)子、重離子與物質(zhì)相互作用的物理機(jī)制復(fù)雜，涉及反應(yīng)眾多，產(chǎn)生了多種次級粒子。在這里主要考慮對劑量影響排在前4位的次級粒子，正如質(zhì)子與物質(zhì)相互作用僅考慮次級質(zhì)子、次級電子、氘核和α粒子所沉積下的劑量，可加快模擬質(zhì)子、重離子與物質(zhì)相互作用過程的速度，且獲得結(jié)果也在臨床可接受的范圍內(nèi)。

另外，本發(fā)明中考慮了磁場對次級帶電粒子的影響，修正了次級帶電粒子的運(yùn)動方向，進(jìn)而獲得較為準(zhǔn)確的磁場下次級帶電粒子的運(yùn)動軌跡，有助于優(yōu)化原先的蒙特卡洛方法，拓寬它的應(yīng)用范圍。本算法可用于核磁實(shí)時(shí)引導(dǎo)質(zhì)子放射治療或核磁實(shí)時(shí)引導(dǎo)重離子的治療計(jì)劃系統(tǒng)的劑量計(jì)算。

步驟5、基于gpu快速原子加法統(tǒng)計(jì)劑量結(jié)果：

步驟5.1、將t個(gè)批次所獲得的放射源粒子和次級粒子的沉積劑量以三維矩陣的形式存入gpu的全局內(nèi)存中，在所述三維矩陣中的任意元素記為dosei；dosei表示在三維空間中第i個(gè)位置上的劑量；

步驟5.2、獲取gpu中同一個(gè)線程包內(nèi)的有效線程，并將所要更新的全局內(nèi)存地址相同的有效線程從所述同一線程包中篩選出來；對于第k個(gè)有效線程，篩選出同一個(gè)線程包內(nèi)與第k個(gè)線程所要更新的全局內(nèi)存地址相同的其他所有有效線程。

有效線程(activethread)指的是對于程序的分支結(jié)構(gòu)正處于當(dāng)前分支內(nèi)的線程。比如對于一塊if-else代碼，當(dāng)程序正在執(zhí)行else時(shí)，線程包里的32個(gè)線程有些可能本應(yīng)執(zhí)行if分支，這些線程會被硬件屏蔽掉，成為無效線程，剩下就是有效線程。

對于第k個(gè)有效線程，用固有函數(shù)__ballot()篩選出其他所有有效線程，得到一個(gè)32比特的標(biāo)記mk，第i個(gè)比特表示第i個(gè)線程是否是有效線程，1為有效，0為無效。

將第j個(gè)有效線程所要更新的全局內(nèi)存地址aj用固有函數(shù)__shuffle()傳到第k個(gè)有效線程，再用__ballot()篩選出線程包內(nèi)滿足ak＝aj的所有有效線程，將結(jié)果存儲在標(biāo)記mk中：如果ak＝aj，則令mk中第j個(gè)比特為1，否則令其為0。如果所有有效線程均不滿足ak＝aj，則對于第j+1個(gè)有效線程重復(fù)本步操作，直到所有有效線程都處理完畢。

步驟5.3、將篩選出來的有效線程所對應(yīng)的全局內(nèi)存中所存儲的劑量進(jìn)行累加，得到的劑量結(jié)果存入序號最小的有效線程中；

步驟5.4、利用“卡漢求和”以及“比較-交換”算法將所有屬于同一線程包內(nèi)序號最小的有效線程中的劑量結(jié)果累加到gpu的全局內(nèi)存中，從而得到人體模型中放射源粒子和次級粒子的總劑量；

步驟5.5、將所述總劑量除以放射源粒子的總數(shù)n，從而得到歸一化的劑量結(jié)果，并將其返回到cpu內(nèi)存。

以上步驟，從軟件的角度實(shí)現(xiàn)了基于kepler構(gòu)架和maxwell構(gòu)架的gpu支持雙精度原子加法運(yùn)算，解決了基于kepler構(gòu)架和maxwell構(gòu)架的gpu硬件上不支持雙精度原子加法運(yùn)算的問題。具體的是解決線程包內(nèi)部的線程競爭，實(shí)現(xiàn)高精度與低精度的混合計(jì)算時(shí)不丟失高精度的數(shù)據(jù)信息，進(jìn)而提升單精度原子加法運(yùn)算的精度和雙精度原子加法運(yùn)算的速度。

本算法的整個(gè)流程就是利用gpu硬件平臺，采用改進(jìn)的蒙特卡洛方法快速準(zhǔn)確地模擬計(jì)算出病人進(jìn)行質(zhì)子與重離子放射治療或核磁實(shí)時(shí)引導(dǎo)質(zhì)子與重離子放射治療時(shí)沉積的劑量。