CN110694186B - Computer readable storage medium for fitting particle beam spot position based on hardware gauss - Google Patents
Computer readable storage medium for fitting particle beam spot position based on hardware gauss Download PDFInfo
- Publication number
- CN110694186B CN110694186B CN201910893001.6A CN201910893001A CN110694186B CN 110694186 B CN110694186 B CN 110694186B CN 201910893001 A CN201910893001 A CN 201910893001A CN 110694186 B CN110694186 B CN 110694186B
- Authority
- CN
- China
- Prior art keywords
- strip
- fitting
- beam spot
- signal
- ionization chamber
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
- A61N5/00—Radiation therapy
- A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
- A61N5/1048—Monitoring, verifying, controlling systems and methods
- A61N5/1064—Monitoring, verifying, controlling systems and methods for adjusting radiation treatment in response to monitoring
- A61N5/1065—Beam adjustment
Abstract
The invention discloses a computer-readable storage medium for fitting the beam spot position of a particle beam based on hardware gauss, and belongs to the technical field of beam detection. The method comprises the following steps: determining the position range of a fitting strip according to a current amplitude signal sampled on a cathode strip of an ionization chamber; filling the missing values of the strip signals in the fitting strip position range; and fitting by using a Gaussian fitting algorithm in hardware according to the strip signals and the positions after the default filling processing to obtain the positions of the beam spots in the particle therapy. According to the invention, by researching the collection condition of strip charges in the ionization chamber, a Gaussian fitting algorithm of limited range, threshold setting and linear interpolation is applied to the calculation of the beam spot position in treatment and is realized in FPGA hardware, so that the requirement of rapidly and accurately calculating the beam spot position in real time is realized. And a student t distribution function is introduced, the beam spot position is obtained by adopting a maximum likelihood method, and the beam spot position is verified with the FPGA calculation result, so that the accuracy of beam spot position monitoring is further ensured.
Description
Technical Field
The invention belongs to the technical field of beam detection, and particularly relates to a computer-readable storage medium for fitting the beam spot position of a particle beam based on hardware Gaussian.
Background
The ionization chamber is arranged at the tail end of the treatment head in the active particle treatment scanning process to monitor the beam position irradiated in the treatment process in real time, so that the correct beam dose is irradiated at the correct position. For this application, strip ionization chambers are more suitable than pixel ionization chambers because fewer channels need to be read and the detection speed is faster (Braccii, S., et al., Segmented ionization chambers for beam monitoring in sunlight. model Physics Letters A, 2015.30(17): p.1540026.). Since the uniformity of the delivered dose distribution is directly dependent on the beam spot lateral position accuracy, positional deviations in excess of one millimeter can result in dose fluctuations of several percent, due toThis required position accuracy must be on the sub-millimeter scale (actions, o., d.meer, and S).
Journal of Instrumentation 2014.9(12): p.C12037-C12037.) is used. It is therefore important to verify the position of the beam spot quickly and accurately, otherwise the beam may hit the wrong position or the beam may not be turned off in time.
As shown in fig. 1, when a beam vertically strikes the ionization chamber, the beam envelope is gaussian distributed, and the beam passing through the ionization chamber collides with the gas in the ionization chamber, which causes gas molecules to be ionized to generate electron-ion pairs, and the negative high voltage applied to the cathode strip of the ionization chamber attracts positive ions to rapidly separate the electron-ion pairs, thereby inducing a current amplitude signal on the cathode strip. The more the gas molecules are ionized at the position closer to the beam, the stronger the current amplitude signal induced on the strip at the corresponding position is, and thus the amplitude of the current amplitude signal on the strip is also in Gaussian distribution on the X-Y plane. The current amplitude signals on the ionization chamber strips are processed to obtain beam information such as the position of a beam spot, the size of the beam spot and the like. The current calculation method for the beam spot position in the ionization chamber mainly comprises a centroid method and a Gaussian fitting method.
Although the centroid method is simple to calculate, the following drawbacks exist: under the condition of poor signal-to-noise ratio, the fitting accuracy cannot meet the requirement; only the position information of the beam spot can be obtained, and other information such as the size of the beam spot cannot be obtained; in the case of strip damage to the ionization chamber, the centroid method fails; when the beam strikes the edge of the ionization chamber, the complete beam spot information cannot be collected by the strip on the ionization chamber, and the centroid method fails. The traditional Gaussian fitting method has the following defects: under the condition of poor signal-to-noise ratio, the fitting accuracy cannot meet the requirement; in the case of a strip damage in the ionization chamber, the gaussian fitting algorithm fails; when the beam is shot at the edge of the ionization chamber, the complete beam spot information cannot be collected by the strip on the ionization chamber, and the Gaussian fitting method is invalid; when the method is realized in the FPGA, the calculated amount is large, and the real-time monitoring of the beam spot position is not facilitated.
Disclosure of Invention
In view of the shortcomings and needs in the art, the present invention provides a computer-readable storage medium for hardware-based Gaussian fitting of beam spot position of a particle beam, which aims to quickly and accurately fit the beam spot position information on a strip ionization chamber.
To achieve the above object, according to a first aspect of the present invention, there is provided a beam spot position method of a particle beam therapy device based on hardware gaussian fitting, the method comprising the steps of:
s1, determining a position range of a fitting strip according to a current amplitude signal sampled on a cathode strip of an ionization chamber;
s2, filling strip signal missing values in the fitting strip position range;
s3, fitting the strip signals and the positions after the default filling processing in hardware by using a Gaussian fitting algorithm to obtain the position lambda of the beam spot in the particle therapy1。
Specifically, step S1 includes the following sub-steps:
s11, according to the amplitude x of each strip current amplitude signal sampled on the cathode strip of the ionization chamberiFinding the peak x of the band signalmax；
S12, according to the peak value x of the strip signalmaxSetting a strip signal threshold Thres ═ x for position fittingmaxX alpha, the first preset value alpha epsilon (0, 1);
s13, judging the current amplitude signal amplitude x of each stripiIf it is higher than the threshold Thres, determining the fitted band range [ i [ ]min,imax]。
Specifically, step S2 includes the following sub-steps:
s21, according to the strip signal peak value xmaxSetting a threshold value D ═ x for linear interpolationmaxX beta, and a second preset value beta epsilon (0, 1);
s22, judging whether to use a threshold D of linear interpolation
i∈[imin,imax]If yes, for the signal value xiLinear interpolation is performed.
Specifically, step S3 includes the following sub-steps:
s31, calculating in an FPGA to obtain a first matrix according to the strip position i after the default value filling processing
Wherein N is the number of the sample points after pretreatment, and i is the position of the ionization chamber strip;
s32, according to the strip signal x and the position i after the missing value filling processing, a CORDIC method is adopted in the FPGA to obtain a second matrix through calculation
Wherein x is the signal amplitude on the cathode strip;
s33, according to the first matrix and the second matrix, matrix triangular decomposition is adopted in the FPGA, and a matrix equation containing unknown parameters a, b and c is solved;
s34, solving the position of the beam spot according to the values of the parameters c and b
Specifically, the initial value of CORDIC iteration is set to x0＝(w+1)·2m，y0＝(w-1)·2m，z0And (5) converting the floating-point number operation into an integer operation when the floating-point number operation is 0.
Specifically, the method further comprises the following steps:
s4, setting data with the current amplitude signal smaller than 0 collected on the ionization chamber strip to zero;
s5, constructing a strip signal sample corresponding to the number frequency of positive ions according to the strip current signal amplitude after the zero setting treatment;
s6, constructing an objective function based on the student t process of generating the sample by the strip signal;
s7, carrying out parameter estimation on the target function through a maximum likelihood method to obtain the calculated beam spot central position lambda2；
S8, comparing lambda1And λ2And checking the beam spot position.
Specifically, for the signal value { x1,x2,x3,…,xnRounding off to construct band signal samples corresponding to the number frequency of positive ions
Specifically, the constructed objective function f (i | λ)2ρ, ν) is as follows:
wherein i is the position of the strip, λ2And (4) the central position of the beam spot, rho is the radius of the beam spot, upsilon is the degree of freedom of distribution, and tau (·) is a gamma function.
Specifically, in step S8, if | λ1-λ2If the | is less than or equal to k, the accuracy of the result of fitting by using a Gaussian fitting algorithm in hardware is considered to be within an acceptable range, and the treatment is continued; otherwise, the fitting result is deemed to be unreliable, the treatment is discontinued, the check is performed again, and k is the third preset value.
To achieve the above object, according to a second aspect of the present invention, there is provided a computer-readable storage medium having stored thereon a computer program which, when executed by a processor, implements the beam spot position method of a hardware-based gaussian-fit particle beam therapy apparatus according to the first aspect.
Generally, by the above technical solution conceived by the present invention, the following beneficial effects can be obtained:
(1) according to the invention, by researching the collection condition that the beam is irradiated on the ionization chamber and the strip in the ionization chamber has charges, an improved Gaussian fitting algorithm of limited range, threshold setting and linear interpolation is applied to the calculation of the beam spot position in treatment, and is realized in FPGA hardware, so that the requirement of quickly and accurately calculating the beam spot position in real time is realized.
(2) The invention carries out data preprocessing on the current amplitude signal on the ionization chamber strip to generate a strip signal sample corresponding to the number frequency of positive ions, adopts a maximum likelihood method to obtain the beam spot position after t distribution is introduced according to the sample, and checks the beam spot position with the FPGA calculation result, thereby further ensuring the accuracy of beam spot position monitoring.
Drawings
FIG. 1 is a schematic diagram of the operation of a prior art strip ionization chamber;
FIG. 2 is a flowchart of a beam spot positioning method for a hardware-based Gaussian fit particle beam therapy device according to an embodiment of the present invention;
FIG. 3 is a graph of a Gaussian fit result without a fitting threshold provided by an embodiment of the present invention;
FIG. 4 is a graph of a Gaussian fit result with a 20% fitting threshold set according to an embodiment of the present invention;
fig. 5 is a graph of fitting results of t-distribution of students provided by the embodiment of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention. In addition, the technical features involved in the embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.
The invention conception of the invention is as follows: the traditional Gaussian fitting algorithm is improved and is realized in an FPGA. The collected induced current amplitude signals on the ionization chamber strips are preprocessed to form samples, the samples are modeled by introducing a student t distribution function, compared with a Gaussian function, the student t distribution function has longer trailing, the influence of terminal strip signals (namely outliers) seriously influenced by noise can be effectively inhibited, and the fitting accuracy is high under the condition of small samples. The fitting efficiency and the accuracy are greatly improved, the beam spot position is calculated by using a maximum likelihood method, and the online and offline verification is carried out, so that the treatment safety is further ensured. The requirement of real-time on-line beam position monitoring is met, and submillimeter-level precision is realized.
As shown in fig. 2, the present invention provides a spot location method of a particle beam therapy device based on hardware gaussian fitting, the method comprising the steps of:
s1, determining the position range of a fitting strip according to a current amplitude signal sampled on a cathode strip of an ionization chamber.
And (3) according to the collected cathode strip signal of the ionization chamber, adopting an improved Gaussian fitting algorithm in hardware to monitor the position of the beam spot in real time. If the hardware onboard clock frequency is higher than 1.5MHz, for example, the DSP is suitable for a system with a lower sampling rate, a low data rate, multi-conditional operation, complex multi-algorithm task processing, C language programming and floating point, and is suitable for multi-conditional processing at a lower sampling rate, particularly complex multi-algorithm tasks. The FPAG is suitable for high-speed sampling rate (MHZ level) of a system, high data rate, block diagram mode programming, fixed or repeated processing tasks and fixed point use, and is suitable for the situation of high-speed sampling frequency, particularly the situation of fixed or repeated tasks, and the situation of trial-manufacture prototype and system development. Since the FPGA is suitable for the high-speed sampling rate (MHZ level) of the system, the FPGA is selected in the embodiment, and the calculation time can be less than 100 μ s when the on-board clock frequency is higher than 1.5 MHZ.
Under the condition of considering FPGA computing resources and computing precision, converting the acquired strip current amplitude signal sample data into 16-bit binary data, wherein 8 bits are integer bits, and the rest are decimal bits.
Because the size of the beam spot in the active scanning particle therapy is only a few millimeters, and the width of a strip of a beam monitor-strip ionization chamber which is vertically arranged with the beam is generally 1-3 millimeters, when the beam passes through the ionization chamber to ionize gas, only a few cathode strips near the beam can sense effective current amplitude signals, and the current amplitude signals acquired by the rest strips far away from the center of the beam spot are weak and are seriously interfered by noise. Based on the situation, preprocessing is carried out on the collected current amplitude signal amplitude samples on the cathode strip of the ionization chamber, the position of a peak value is determined, a fitting threshold value is set, and a weak signal strip far away from the beam spot position is filtered. The tail strip signals seriously interfered by the substrate noise can be filtered, so that the influence of the noise on the fitting signals is inhibited, and the signal-to-noise ratio is greatly improved.
S11, according to the amplitude x of each strip current amplitude signal sampled on the cathode strip of the ionization chamberiFinding the peak x of the band signalmax。
The maximum current amplitude signal x collected in the strip ionization chamber for dose monitoring can be found by adopting a divide-and-conquer methodmaxA corresponding strip.
S12, according to the peak value x of the strip signalmaxSetting a strip signal threshold Thres ═ x for position fittingmaxX α, the first preset value α ∈ (0, 1).
By fitting a threshold value Thres, weak strip current amplitude signals away from the beam spot location can be filtered out. Because the fitting range is reduced to the left and the right strips of the central position of the beam spot, the number of fitted strips is reduced, the calculated amount is reduced, and the calculation efficiency is improved. In the present embodiment, the first preset value α is 20%.
S13, judging the current amplitude signal amplitude x of each stripiIf it is higher than the threshold Thres, determining the fitted band range [ i [ ]min,imax]。
Judging the current amplitude signal amplitude x of each stripiIf it is higher than the threshold Thres, and if so, its state value sign [ i ] is marked]Is 1, otherwise, its state value sign [ i ] is marked]Is 0. At the peak x of the band signalmaxIn the left and right range of the stripe, find out the satisfied state value sign [ i ]]Leftmost strip i of 1minAnd the rightmost strip imax。
And S2, filling the strip signal missing values in the fitting strip position range.
Aiming at the situation that partial strips possibly exist in the actual situation are damaged and cannot work normally, so that corresponding strip signals are lost, the method judges whether a damaged cathode strip exists in the range of the irradiated beam spot, namely the strip has no induced current amplitude signal or the signal amplitude is extremely weak and is far lower than the signals on the strips at two sides, and if the situation exists, the strip is filled, so that the accuracy of a calculation result is ensured.
S21, according to the strip signal peak value xmaxSetting a threshold value D ═ x for linear interpolationmaxX β, the second preset value β ∈ (0, 1).
The linear interpolation is simple and convenient, and hardware implementation is facilitated, so the linear interpolation is preferred by the filling method of the embodiment. And judging whether an abnormal band with a low signal amplitude exists in the fitting range through a threshold value D of linear interpolation, wherein the second preset value beta is 5% in the embodiment.
S22, judging whether to use a threshold D of linear interpolation
i∈[imin,imax]If yes, for the signal value xiLinear interpolation is performed.
xi<D, indicating that there is a band damage condition near the beam spot position, for signal value xiLinear interpolation is performed, and the interpolation formula is as follows:
s3, fitting the strip signals and the positions after the default filling processing in hardware by using a Gaussian fitting algorithm to obtain the position lambda of the beam spot in the particle therapy1。
After sample data of all strip current signal amplitudes are processed, the sample data, namely strip signal values, of a plurality of effective strip current amplitude signals only containing the central position of a beam spot is obtained
And its position { imin,…,imaxAnd performing Gaussian fitting calculation on the data to obtain the position of the beam spot in the particle therapy.
Since the beam spot is gaussian distributed on the X, Y plane, the induced signal values of the X, Y direction stripe should also be gaussian distributed, and the stripe signal value x and the stripe position i satisfy the gaussian distribution.
For a 1-dimensional random variable i, the gaussian distribution function is given by:
wherein, x is the amplitude of the strip current signal, A is the peak value of the strip signal, lambda is the central position of the beam spot, rho is the radius of the beam spot, and i is the position of the strip.
The above-described non-linearized equation with unknown parameters is changed to a linear equation. Taking logarithm of the above equation to obtain
Wherein the content of the first and second substances,
the error function is
(i)＝lnx-(a+bi+ci2)
Wherein, (i) is the fitting value obtained by the amplitude x of the signal on each strip and under different parameters a, b and c
The difference between the values of the two, wherein,
using least square method to process the strip signal sample-strip signal value
And its position { imin,…,imaxAnd performing Gaussian fitting to obtain the position of the beam spot. The objective is to find the values of the parameters a, b, c that minimize the value of the error function (i).
S31, calculating in an FPGA to obtain a first matrix according to the strip position i after the default value filling processing
Wherein N is the number of the sample points after pretreatment, and i is the position of the ionization chamber strip.
S32, according to the strip signal x and the position i after the missing value filling processing, a CORDIC method is adopted in the FPGA to obtain a second matrix through calculation
Where x is the signal amplitude on the cathode strip.
Since logarithmic calculation cannot be directly performed in the FPGA, the objective function is approximated by shift and addition operations using CORDIC (coordinate rotation digital calculation) calculation technology. The CORDIC function in hyperbolic coordinate system is as follows:
xn+1＝xn+dn(2-nyn)
yn+1＝yn+dn(2-nxn)
zn+1＝zn-dntanh-12-n
wherein d isn＝-sign(xnyn)。
In addition, convergence can be satisfied only when a specific iteration of n — 4,13, …, k,3k +1, …, etc. is repeated. For FPGA, calculate 2nIs advantageous because it only needs to move n bits left or right according to the sign of n, and tanh can be stored in the FPGA in advance -12-nThe value of (c).
Taking an initial iteration value x0＝w+1，y0＝w-1，z0When 0, after the iteration is completed, the following equation can be obtained:
to solve the problem of CORDIC function convergence domain, the logarithm of the band signal value x is decomposed as follows:
lnx＝ln(w·2E)＝Eln2+lnw
in order to save computing resources while ensuring computing precision, an initial value of CORDIC iteration is set to be x0＝(w+1)·2m，y0＝(w-1)·2m，z0And when the floating-point number operation is converted into the integer operation, the floating-point number operation is successfully avoided, and the calculation resource in the FPGA is saved while the calculation precision is ensured. Is calculated to
Substituting the above equation, the logarithmic value lnx of the band signal is finally obtained.
For example, x is 10, ln10 is ln (1.25 · 2)3)＝3ln2+ln1.25
Setting an initial value x of an iteration0＝(1.25+1)·210＝2304，y0＝(1.25-1)·210＝256，z 00. To satisfy convergence, certain iterations of n-4, 13, …, k,3k +1, …, etc. are repeated. Z is obtained through 16 times of iterative calculation160.1116, ln1.25 ═ 2 × z16＝0.2231，ln10＝ln(1.25·23)＝3ln2+ln1.25＝3·0.6931+0.2231＝2.3024。
S33, according to the first matrix and the second matrix, matrix triangular decomposition is adopted in the FPGA to solve a matrix equation containing unknown parameters
To calculate the parameters a, b, c in the system of linear equations. Requiring taking the matrix AInverse matrix, so an equation like AX ═ B can be solved as X ═ a-1Form B.
By triangulating the matrix, the equation can be simplified as follows:
wherein A ═ P-1U is LU, L is a lower triangular matrix, and U is an upper triangular matrix. Thus, the equation AX ═ LUX ═ B can be solved as:
UX＝PB
s34, solving the position lambda of the beam spot according to the values of the parameters c and b1。
Since the measured strip current amplitude signals are spatially gaussian distributed, by fitting the strip current amplitude signals in the X, Y directions separately, ix＝λx，iy＝λyObtaining the central position (i) of the beam spotx，iy). The two-dimensional gaussian function problem is simplified to a one-dimensional gaussian function problem. Specifically, a least square method is used for carrying out iterative calculation on unknown parameters in the Gaussian function, and finally the central position of the beam spot is obtained. The method comprises the following steps of carrying out arithmetic such as logarithm taking and matrix inversion on the FPGA.
When a charged particle moves in a medium, it loses energy by interacting with molecules in the medium. The ionization chamber consists essentially of two electrodes with a volume of gas between them, which are typically connected to a high voltage power supply of 100v to 1000v, allowing the chamber to operate in the saturation region. As the proton beam passes through the ionization chamber, the proton beam loses energy by colliding with the gas in the ionization chamber, creating electron-ion pairs. In the absence of an applied drift field, electron-ion pairs rapidly recombine due to recombination effects. The high voltage applied to the plate separates the electrons and ions from each other and the signal produced by the moving charge provides information about the energy and position of the proton beam. In clinical treatment, the ionization chamber is placed perpendicular to the direction of the beam current, collecting charges proportional to the energy deposited in the gas of the ionization chamber.
The strip signal containing certain substrate noise can be obtained through the strip ionization chamber, and the approximate position of the beam spot can be preliminarily determined by finding the peak value of the strip signal, because when the beam passes through the ionization chamber, the ionization degree of protons and gas molecules in the ionization chamber is highest on the path through which the beam passes, the charge collected by the adjacent strip is the largest, and the generated current amplitude signal is strongest. However, merely determining which band the spot is centered on is far from accurate for patient treatment, and further processing of the band signals is required to achieve sub-millimeter accuracy resolution, from safety considerations.
For active scanning, in order to realize conformal treatment, cancer cells are killed to the maximum extent and normal tissues are not damaged as much as possible, the beam current is weak, the magnitude is nA, only a few bands can collect charges after ionization, and the end band signals far away from the beam are extremely weak and are seriously influenced by substrate noise. By analyzing the actual situation, the fitting range is limited by taking the peak value strip as the center, so that the substrate noise is suppressed, and the signal-to-noise ratio is improved. The fitting results and their alignment are shown in fig. 3 and 4. Fig. 3 shows a gaussian fitting result without a fitting threshold, which shows that the fitting curve is seriously interfered by tail noise and the fitting result is poor; fig. 4 is a gaussian fitting result with a fitting threshold of 20%, and it can be seen that the fitting result is obviously improved after tail noise is filtered.
Meanwhile, aiming at the situation that a certain strip in the ionization chamber can not work normally and the signal of the strip is lost, which can happen in the real situation, the signal is in a defined fitting range (i)min,imax) And linear interpolation is carried out on unreliable strip signals, so that the calculation accuracy is further improved.
Example one
The method further comprises the following steps:
and S4, setting the data with the current amplitude signal smaller than 0 collected on the ionization chamber strip to zero.
For the ionization chamber strip {1,2,3, …, n }, the corresponding induced current amplitude signal value is { x }1,x2,x3,…,xnV. current amplitude signal value xiWhether to judge in proportion to the number of positive ions attracted by the cathode strip i, i.e., the frequency of migration of the positive ions to the strip i
i∈[1，n]If yes, let xi＝0。xi<0 indicates that the signal is mainly substrate noise, and the number of positive ions attracted by the strip is little or almost none, so that the fitting error caused by the substrate noise is reduced on one hand, and the subsequent data processing is facilitated on the other hand.
And S5, constructing a strip signal sample corresponding to the number frequency of the positive ions according to the strip current signal amplitude after the zero setting treatment.
Based on the operating principle of the ionization chamber, the higher the amplitude of the strip current signal, the more particles incident at the position, the more positive ions collected on the strip, that is, the more positive ions are attracted by the electric field and migrate to the strip, and from the viewpoint of probability theory, the amplitude of the strip current signal is proportional to the frequency of migration of the positive ions to the strip.
For the signal value { x1,x2,x3,…,xnRound off to generate samples
For example, if corresponding to strips 1,2,3 … n, there is a signal value x1＝0.4，x2＝0.8，x3＝1.2，…，xnAt 7.3, the signal value is rounded by {0,1,1, …,7}, and a sample is generated
This indicates that the number of positive ions reaching the band position 2 in the generated sample is 1, and the number of positive ions reaching the band position 3 is 1, and thusAnd so on.
S6, constructing an objective function f (i | lambda, rho, upsilon) by a student t process based on the strip signal generation sample:
compared with a Gaussian function, the student t distribution function has longer tailing, can effectively inhibit the influence of terminal strip signals (namely outliers) seriously influenced by noise, and has high fitting accuracy under the condition of small samples.
For a 1-dimensional random variable i, the probability density function of the Student-t distribution is given by:
where μ is the sample mean, σ is the sample standard deviation, and upsilon is the degree of freedom of the distribution, and when the degree of freedom upsilon → ∞, the student t distribution converges to a gaussian distribution having the same mean and variance. τ (i) is the gamma function, in the real domain:
s7, carrying out parameter estimation on the target function through a maximum likelihood method to obtain the calculated beam spot central position lambda2。
For generating samples
Capacity of it
The combined probability density function is
Constructing likelihood functions of samples
Logarithm of likelihood function
And solving partial derivatives of the equation, enabling the equation to be equal to 0, and solving the maximum likelihood estimation of the beam spot position lambda.
S8, comparing lambda1And λ2And checking the beam spot position.
If lambda1-λ2If the | is less than or equal to k, the accuracy of the result of fitting by using a Gaussian fitting algorithm in hardware is considered to be within an acceptable range, and the treatment is continued; otherwise, the fit is deemed to be unreliable, treatment is discontinued, and re-calibration is performed. The value range of k in this example is [0,0.5mm ]]。
An mle (maximum likelihood estimation) function can be used in MATLAB to directly solve unknown parameters in t distribution to obtain the central position lambda of the beam spot2。
Data samples
Due to the constructed sample
Is from a normal distribution N (mu, sigma)2) A sample of
Wherein the content of the first and second substances,
and s2Sample mean and sample variance, respectively.
Mle (a, 'distribution', 'tlocationscale'); % returns the maximum likelihood estimate, phat, for the specified distribution, tlocatationscale.
Obtaining the central position lambda of the beam spot2α t (1), α t (2) is the beam spot radius, and α t (3) is the degree of freedom.
The sample data is preprocessed, and a fitted t distribution curve is shown in fig. 5, so that the tail interference signal is effectively filtered due to the characteristic of t distribution and heavy tail.
It will be understood by those skilled in the art that the foregoing is only a preferred embodiment of the present invention, and is not intended to limit the invention, and that any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the scope of the present invention.
Claims (6)
1. Computer-readable storage medium for hardware-gaussian-based fitting of a beam spot position of a particle beam, characterized in that the computer-readable storage medium has stored thereon a computer program which, when executed by a processor, implements a method for hardware-gaussian-based fitting of a beam spot position of a particle beam therapy device, the method comprising the steps of:
s1, determining a position range of a fitting strip according to a current amplitude signal sampled on a cathode strip of an ionization chamber;
s2, filling strip signal missing values in the fitting strip position range;
s3, fitting the strip signals and the positions after the default filling processing in hardware by using a Gaussian fitting algorithm to obtain the position lambda of the beam spot in the particle therapy1；
Step S1 includes the following substeps:
s11, according to the amplitude x of each strip current amplitude signal sampled on the cathode strip of the ionization chamberiFinding the peak x of the band signalmax；
S12, according to the peak value x of the strip signalmaxSetting a strip signal threshold Thres ═ x for position fittingmaxX alpha, the first preset value alpha epsilon (0, 1);
s13, judging the current amplitude signal amplitude x of each stripiIf it is higher than the threshold Thres, determining the fitted band range [ i [ ]min，imax]；
Step S2 includes the following substeps:
s21, according to the strip signal peak value xmaxSetting a threshold value D ═ x for linear interpolationmaxX beta, and a second preset value beta epsilon (0, 1);
s22, judging whether to use a threshold D of linear interpolation
If so, for the signal value xiLinear interpolation is performed.
2. The computer-readable storage medium of claim 1, wherein step S3 includes the sub-steps of:
s31, calculating in an FPGA to obtain a first matrix according to the strip position i after the default value filling processing
Wherein N is the number of the sample points after pretreatment, and i is the position of the ionization chamber strip;
s32, according to the strip signals and the position i after the missing value filling processing, a CORDIC method is adopted in the FPGA to obtain a second matrix through calculation
Wherein x is the signal amplitude on the cathode strip;
s33, according to the first matrix and the second matrix, matrix triangular decomposition is adopted in the FPGA, and a matrix equation containing unknown parameters a, b and c is solved;
s34, according to the parameters c andb, solving the position of the beam spot
3. The computer-readable storage medium of claim 1 or 2, wherein the method further comprises the steps of:
s4, setting data with the current amplitude signal smaller than 0 collected on the ionization chamber strip to zero;
s5, constructing a strip signal sample corresponding to the number frequency of positive ions according to the strip current signal amplitude after the zero setting treatment;
s6, constructing an objective function based on the student t process of generating the sample by the strip signal;
s7, carrying out parameter estimation on the target function through a maximum likelihood method to obtain the calculated beam spot central position lambda2；
S8, comparing lambda1And λ2And checking the beam spot position.
4. The computer-readable storage medium of claim 3, wherein { x } is for a signal value1，x2，x3，...，xnRounding off to construct band signal samples corresponding to the number frequency of positive ions
5. The computer-readable storage medium of claim 3, wherein the constructed objective function f (i | λ)2ρ, ν) is as follows:
where i is the ionization chamber strip position, λ2Is the central position of the beam spot, rho is the radius of the beam spot, upsilon is the degree of freedom of distribution, tau (·)Is a gamma function.
6. The computer-readable storage medium of claim 4 or 5, wherein in step S8, if λ1-λ2If the | is less than or equal to k, the accuracy of the result of fitting by using a Gaussian fitting algorithm in hardware is considered to be within an acceptable range, and the treatment is continued; otherwise, the fitting result is deemed to be unreliable, the treatment is discontinued, the check is performed again, and k is the third preset value.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910893001.6A CN110694186B (en) | 2019-09-20 | 2019-09-20 | Computer readable storage medium for fitting particle beam spot position based on hardware gauss |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910893001.6A CN110694186B (en) | 2019-09-20 | 2019-09-20 | Computer readable storage medium for fitting particle beam spot position based on hardware gauss |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110694186A CN110694186A (en) | 2020-01-17 |
| CN110694186B true CN110694186B (en) | 2020-11-17 |
Family
ID=69196271
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910893001.6A Active CN110694186B (en) | 2019-09-20 | 2019-09-20 | Computer readable storage medium for fitting particle beam spot position based on hardware gauss |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110694186B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112698380A (en) * | 2020-12-16 | 2021-04-23 | 南京大学 | Beam section processing method suitable for low-energy proton beam under strong background noise |
| CN115040792B (en) * | 2022-03-25 | 2023-03-07 | 中国原子能科学研究院 | Signal generating device of ionization chamber for proton treatment |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP1852717A1 (en) * | 2006-02-10 | 2007-11-07 | Commissariat A L'energie Atomique | Method of estimating scattered radiation in a two-dimensional detector |
| CN102858407A (en) * | 2010-02-10 | 2013-01-02 | 株式会社东芝 | Particle Beam Irradiation Device And Control Method Therefor |
| CN104338244A (en) * | 2013-08-07 | 2015-02-11 | 株式会社日立制作所 | Beam monitor system and particle beam irradiation system |
| CN105824817A (en) * | 2015-01-05 | 2016-08-03 | 苏州瑞派宁科技有限公司 | Flash pulse digitization method |
| CN107462918A (en) * | 2017-08-22 | 2017-12-12 | 合肥中科离子医学技术装备有限公司 | A kind of accelerator beam cross section measuring system and method based on LabVIEW |
| CN108681992A (en) * | 2018-04-23 | 2018-10-19 | 南京理工大学 | The image interpolation algorithm of laser facula is measured for detector array method |
| US10191160B1 (en) * | 2018-08-31 | 2019-01-29 | David Edward Newman | Staggered detector array for locating radioactive sources |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB0014096D0 (en) * | 2000-06-09 | 2000-08-02 | Council Cent Lab Res Councils | Proportional gas counters |
| JP4151703B2 (en) * | 2006-04-04 | 2008-09-17 | 日新イオン機器株式会社 | Ion beam measuring apparatus, measuring method, and ion beam irradiation apparatus |
| JP5668000B2 (en) * | 2012-03-02 | 2015-02-12 | 株式会社日立製作所 | Beam monitor system and particle beam irradiation system |
| CN108671418A (en) * | 2018-05-24 | 2018-10-19 | 中国科学院近代物理研究所 | Guide of magnetic resonant image device for ion beam radiation therapy |
-
2019
- 2019-09-20 CN CN201910893001.6A patent/CN110694186B/en active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP1852717A1 (en) * | 2006-02-10 | 2007-11-07 | Commissariat A L'energie Atomique | Method of estimating scattered radiation in a two-dimensional detector |
| CN102858407A (en) * | 2010-02-10 | 2013-01-02 | 株式会社东芝 | Particle Beam Irradiation Device And Control Method Therefor |
| CN104338244A (en) * | 2013-08-07 | 2015-02-11 | 株式会社日立制作所 | Beam monitor system and particle beam irradiation system |
| CN105824817A (en) * | 2015-01-05 | 2016-08-03 | 苏州瑞派宁科技有限公司 | Flash pulse digitization method |
| CN107462918A (en) * | 2017-08-22 | 2017-12-12 | 合肥中科离子医学技术装备有限公司 | A kind of accelerator beam cross section measuring system and method based on LabVIEW |
| CN108681992A (en) * | 2018-04-23 | 2018-10-19 | 南京理工大学 | The image interpolation algorithm of laser facula is measured for detector array method |
| US10191160B1 (en) * | 2018-08-31 | 2019-01-29 | David Edward Newman | Staggered detector array for locating radioactive sources |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110694186A (en) | 2020-01-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Bassompierre et al. | First determination of the proton electromagnetic form factors at the threshold of the time-like region | |
| EP1289627B1 (en) | Nanodosimeter based on single ion detection | |
| Brown et al. | Coincidence electroproduction of charged pions and the pion form factor | |
| CN110694186B (en) | Computer readable storage medium for fitting particle beam spot position based on hardware gauss | |
| Baller | Liquid argon TPC signal formation, signal processing and reconstruction techniques | |
| Pennazio et al. | Carbon ions beam therapy monitoring with the INSIDE in-beam PET | |
| AU2001274814A2 (en) | Nanodosimeter based on single ion detection | |
| Pennazio et al. | Proton therapy monitoring: spatiotemporal emission reconstruction with prompt gamma timing and implementation with PET detectors | |
| Nojiri et al. | Reconstructing particle masses from pairs of decay chains | |
| Va'Vra | High resolution drift chambers | |
| Åkesson et al. | The ATLAS TRT straw proportional tubes: Performance at very high counting rate | |
| Hilgers et al. | Reducing the background of secondary ions in an ion-counting nanodosimeter | |
| Blevins et al. | π+− p Interactions in the Energy Range around 500 Mev | |
| Anderson et al. | Proton total reaction cross section measurements for Ca 4 0, 4 4, 4 8 at 700 MeV | |
| Diwan et al. | A novel method for event reconstruction in Liquid Argon Time Projection Chamber | |
| Blyth et al. | Nonnegative Gaussian process tomography for generalized segmented planar detectors | |
| Lin et al. | Fast and accurate position verification algorithms for dose delivery system in proton therapy scanning nozzle system | |
| Zhang et al. | Simulation and data analysis of a fission time projection chamber based on Monte Carlo method | |
| Pietrini | Statistical processing of Flash X-ray Imaging of protein complexes | |
| Bukharskii et al. | On the proton radiography of magnetic fields in targets irradiated by intense picosecond laser pulses | |
| Casiraghi | Nanodosimetric track structure studies for applications in particle therapy | |
| Bercaw et al. | Elastic Scattering of Negative Pions from O 16 in the Region of the (3, 3) Resonance | |
| Dydak | On distortions of TPC coordinates: inhomogeneities of electric and magnetic field | |
| Płażek et al. | Track finding algorithm using the track road method for the PANDA Forward Tracker | |
| Pizzolotto | Measurement of Lambda polarisation observables at the COSY-TOF spectrometer |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |