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

背景技術(shù)：

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟1.3、獲取質(zhì)子和重離子分別與物質(zhì)發(fā)生反應的核數(shù)據(jù)并進行處理，得到所述核數(shù)據(jù)的宏觀截面數(shù)據(jù)∑，對所述宏觀截面數(shù)據(jù)∑按照能量的高低進行降序排序，得到排序后的質(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ù)和輸運任務(wù)的批次；

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

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

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

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

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

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

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

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

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

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

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

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

步驟3.7、計算第w次抽取的第w個放射源粒子sw在人體模型中進行第u+1次輸運時沉積的劑量

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

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

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

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

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

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

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

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

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

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

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

若為次級¹²c粒子，則按照步驟3.4-步驟3.9進行輸運；

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

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

若為次級¹¹c粒子，則按照步驟3.4-步驟3.9進行輸運；

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

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

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

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

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

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

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

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

附圖說明

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

具體實施方式

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟3.7、計算第w次抽取的第w個放射源粒子sw在人體模型中進行第u+1次輸運時沉積的劑量

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

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

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

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

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

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

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

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

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

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

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

若為次級¹²c粒子，則按照步驟3.4-步驟3.9進行輸運；

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

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

若為次級¹¹c粒子，則按照步驟3.4-步驟3.9進行輸運；

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

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

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

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

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

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

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

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

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

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

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

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

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