JPWO2013129258A1 - Particle beam direction determination system, particle beam direction determination method, and particle beam direction determination computer program - Google Patents

Abstract

A particle beam beam direction determination system is provided that determines a beam direction that is robust to patient setup errors. The beam direction image creation unit 1 creates a plurality of beam direction images having a plurality of electron densities viewed from the positions of a plurality of beam sources. The average power spectrum calculation unit 2? Calculates the average power spectrum for each of the plurality of electron density beam direction images. The inclination calculation unit 3A calculates the inclination of the average power spectrum obtained from the beam direction images having a plurality of electron densities. The beam direction determination unit 3B determines the beam direction of the particle beam based on the inclination.

Description

  The present invention relates to a particle beam direction determination system, a particle beam direction determination method, and a particle beam direction determination computer program.

  In particle beam therapy, such as proton therapy or heavy ion therapy (carbon ions), the dose distribution produced by multiple beams is better defined to reduce the dose to normal tissue and increase the dose to the tumor. Control is required. In particle beam therapy, an accurate deep dose distribution around the target tumor can be obtained by using a compensation filter in the system that forms the deep dose. FIG. 24 shows an example of the configuration of an irradiation field forming apparatus for broad beam irradiation. In this apparatus, the particle beam emitted from the particle beam generator 11 passes between the wobbler electromagnets 12 and 13, passes through the scatterer 14, and then scatters, and then enters the ridge filter 15. The ridge filter 15 creates an enlarged Bragg peak in the particle beam, adjusts the range of the particle beam in the body by the range shifter 16, defines the irradiation range by the multi-leaf collimator 17, and places the particle beam where the particle beam stops by the compensation filter 18. The shape is formed and the target (tumor) 19 is irradiated with a particle beam. Conventionally, the compensation filter 18 is created based on the lateral distribution in the beam direction image of the electron density obtained from the three-dimensional electron density image of the computed tomography (CT) image of each patient, and the dose distribution is determined. . However, if the patient moves during treatment and is in a slightly different posture than the original set posture (occurrence of patient setup error) is not taken into account in the treatment plan, The dose distribution formed in can be decreased in the tumor area or increased in the normal tissue area compared to the planned dose distribution. FIG. 25 shows an example of the correspondence relationship between the created compensation filter 18 and the planned target volume PTV that is the target. In FIG. 25, the “water equivalent thickness” on the vertical axis corresponds to the electron density. The compensation filter 18 stops each particle beam (beam) at the end point of the planned target volume PTV. However, when a patient setup error occurs, the positional relationship between the compensation filter 18 and the planned target volume PTV shifts, and it becomes impossible to match the end of the tumor shape with the end of the dose distribution in the beam direction. In particular, as described above, when the electron density beam direction image has an electron density region that changes sharply with a large amplitude change, such a shift causes the tumor to be irradiated with a low dose, and normal Problems that cause fatal over-irradiation of various organs. Such a problem also occurs when the beam irradiation distance is individually controlled based on the electron direction beam direction image without using the compensation filter.

  In order to solve the over-irradiation problem, Japanese Patent Application Laid-Open No. 2011-212395 (Patent Document 1) arranges and arranges a columnar irradiation field on the inner side of the irradiation target, and sets the state as an initial state to the irradiation target. A particle therapy apparatus using a treatment plan of a treatment planning apparatus that adjusts the arrangement of columnar irradiation fields so that the irradiation dose falls within a predetermined range is disclosed.

  Japanese Patent Laid-Open No. 2011-217902 (Patent Document 2) discloses a technique for providing treatment plan information for suppressing the appearance of a high-line region other than a target region in a particle treatment plan apparatus. In this technique, an irradiation spot is arranged on the outline of the irradiation region, and the irradiation spot is determined so that adjacent irradiation spots are equal to or less than a predetermined set value.

JP 2011-212395 A JP 2011-217902 A

  However, conventional treatment planning devices are robust (some patient setups) that prevent patient setup errors from irradiating tumors at low doses and causing fatal over-irradiation of normal organs. There was nothing that determined the particle beam direction (even if there was an error, neither low-dose irradiation nor excessive irradiation occurred).

  It is an object of the present invention to provide a particle beam beam direction determination system, a particle beam beam direction determination method, and a particle beam beam direction determination computer program that determine a beam direction that is robust against a patient setup error.

  The present invention is directed to a particle beam beam direction determination system that determines a beam direction of a particle beam, which is used in a treatment planning apparatus that generates treatment plan information used in a particle beam irradiation system. The particle beam direction determination system of the present invention includes a beam direction image creation unit, a power spectrum calculation unit, and a beam direction determination unit.

  The beam direction image creation unit uses a three-dimensional electron density image created based on a CT value obtained from a treatment plan CT image of an irradiation target obtained by imaging the irradiation target with a CT imaging apparatus. That is, the beam direction image creation unit assumes the installation positions of a plurality of beam sources at predetermined angular intervals around the three-dimensional electron density image, and moves from the installation positions of the respective beam sources toward the planned target volume image. A two-dimensional image obtained by projecting the electron density from the surface of the irradiation target in the three-dimensional electron density image to the end point of the planned target volume on the virtual two-dimensional screen is installed on the beam source. A beam direction image of electron density viewed from the position is created for each installation position of a plurality of beam sources. Note that the two-dimensional screen is a virtual calculation screen for obtaining a beam direction image by calculation, and the partial projection of the three-dimensional electron density image on the two-dimensional screen is also a calculation projection. That is, the projection of the electron density in the beam direction image creation unit can be realized by sampling and adding the electron density on the beam line from the surface of the irradiation target to each end point of the planned target volume. For example, when there is a cavity in the irradiation target, the electron density of the region including the cavity in the beam direction image of the electron density appears to be low, and for the part (muscle, bone, etc.) where there is no cavity in the irradiation target, The electron density of the region in the electron density beam direction image becomes high.

  A three-dimensional electron density image can be created as follows, for example. First, isotropic voxelization of the treatment plan CT image and the plan target volume in the treatment plan CT image is performed. Secondly, the CT value of the treatment plan CT image converted into isotropic voxels is converted into the relative electron density of water to create a converted electron density image. Third, a three-dimensional electron density image is obtained by extracting a body region (non-air region) portion from a three-dimensional converted electron density image. Using the three-dimensional electron density image thus obtained, the range from the surface of the irradiation target to each end point of the planned target volume is calculated, and based on this range, the electron on the beam line is calculated. A beam direction image of the electron density is created by adding the samples obtained by sampling the density.

  The power spectrum calculation unit calculates a power spectrum for each of the beam direction images having a plurality of electron densities viewed from the positions of the plurality of beam sources. In the present invention, an integrated power spectrum and / or an average power spectrum can be used as the power spectrum. The inventor determines whether there is an electron density region part that changes sharply with a large amplitude change in the beam direction image (determining that the dose distribution generated in such a beam direction is vulnerable to patient setup errors. ) Found that the power spectrum can be utilized. That is, the inventor has found that a power spectrum with many high-frequency components is generated in the spatial frequency component from the region of the electron density that changes sharply with a large amplitude change. Therefore, for example, assuming the arrangement position of the beam source around the 3D electron density image at intervals of 5 ° (if the incident direction of the beam is changed around the 3D electron density image at intervals of 5 °). The average power spectrum is calculated for each of the 72 electron density beam direction images. Specifically, the power spectrum calculation unit first substitutes the average value of the electron density values in the image of the beam direction image of the electron beam into the image peripheral part of the beam direction image, and smoothes the edge of the beam direction image. In the smoothing, a Gaussian filter is applied around the contour of the region having the value of the beam direction image to smooth the edge of the region having the value of the beam direction image. By this smoothing, it is prevented that an originally unnecessary high-frequency component is included in the average power spectrum from an edge generated when an electron density beam direction image is created. Next, preprocessing for normalizing the average value of the entire smoothed image to 0 is performed. Specifically, the average value of the entire image is subtracted from the value of the entire image. Then, a power spectrum image is obtained by calculation based on an image in which the average value is normalized to 0 (preprocessing is performed). The power spectrum image is imaged by performing a two-dimensional Fourier transform on the preprocessed image to obtain a power spectrum by calculation. In other words, the preprocessed electron density beam direction image data is Fourier transformed to obtain power for each frequency bandwidth (frequency resolution) and image it to obtain a power spectrum image. Then, the image areas of the power spectrum image are exchanged so that the zero frequency comes to the center of the image. Next, polar coordinate conversion of the power spectrum is performed to convert the power spectrum image into an image represented by coordinates of spatial frequency and angle. Next, an integrated power spectrum at each spatial frequency of the power spectrum image subjected to polar coordinate conversion is calculated. An average power spectrum is obtained from this integrated power spectrum. The average power spectrum is obtained by taking the natural logarithm of the vertical axis value and the horizontal axis value of the figure showing the integrated power spectrum.

  As described above, the inventor has found that a beam direction image that is vulnerable to a setup error has a power spectrum with a high frequency component in the spatial frequency component. Based on this knowledge, the beam direction determining means determines the beam direction of the particle beam based on the amount of the high frequency component in the spatial frequency component of the power spectrum. That is, the beam direction determining means determines that the direction of the beam source capable of obtaining a beam direction image in which the amount of the high-frequency component is equal to or smaller than a predetermined component amount threshold is an appropriate direction as the beam direction of the particle beam. Therefore, according to the present invention, it is possible to easily determine a beam direction robust to a setup error.

  Whether or not the amount of the high frequency component is equal to or less than a predetermined component amount threshold can be determined based on several evaluation indexes. One evaluation index is the absolute value of the slope of the average power spectrum, and the other evaluation index is a zeroth-order moment obtained from the integrated power spectrum.

  When the absolute value of the slope of the average power spectrum is used as an evaluation index, the power spectrum calculation unit is configured by an average power spectrum calculation unit that calculates an average power spectrum for each of the plurality of beam direction images of electron density. The beam direction determining means includes an inclination calculating unit for calculating an absolute value of the average power spectrum inclination obtained for each of the beam direction images having a plurality of electron densities, and the absolute value of the average power spectrum inclination is equal to or greater than a predetermined inclination threshold value. In such a case, it can be configured that the beam direction determination unit determines that the beam direction of the particle beam is determined by determining that the amount of the high-frequency component is equal to or less than the component amount threshold. The inclination calculation unit calculates the absolute value of the inclination of the average power spectrum obtained from the plurality of electron density beam direction images. The slope or slope of the average power spectrum is obtained as the slope of a first order polynomial with the power spectrum for all frequencies up to the Nyquist frequency. This inclination can be obtained as an inclination of this approximate line by obtaining an approximate line from the average power spectrum by the least square method.

  The beam direction determination unit determines the beam direction of the particle beam based on this inclination. When the region of the electron density that changes sharply with a large amplitude change in the beam direction image becomes smaller, the spatial frequency of the change in the electron direction beam direction image becomes lower and the amplitude becomes smaller. Based on the knowledge that the high-frequency component is reduced and the gradient of the average power spectrum, that is, the absolute value of the gradient is increased, the absolute value of the gradient is obtained. When the absolute value of this slope is large, even if there are some patient setup errors, there are few electron density regions that change sharply with large amplitude changes in the beam direction image. The likelihood of occurrence is reduced. As a result, it is possible to determine a beam direction that is robust against setup errors.

  Therefore, the beam direction determination unit can determine that the direction of the beam source corresponding to the beam direction image having an electron density whose tilt absolute value is larger than the tilt threshold is an appropriate direction as the beam direction of the particle beam. A plurality of types of beam directions are determined as preferred directions depending on how the tilt threshold is determined.

  The beam direction determining unit can determine that the direction of the beam source corresponding to the beam direction image having an electron density whose absolute value of the inclination is smaller than the threshold is an inappropriate direction as the particle beam direction. If it does in this way, what is necessary is just to employ | adopt the direction which a doctor considers preferable from directions other than an inappropriate direction.

  When the 0th-order moment is used as the evaluation index, the power spectrum calculation unit calculates the integrated power spectrum for each of the plurality of electron density beam direction images. The beam direction determining means includes a zero-order moment calculation unit for calculating the zero-order moment from the integral power spectrum obtained for each of the beam direction images having a plurality of electron densities, and the zero-order moment is less than a predetermined moment threshold value. Includes a beam direction determination unit that determines the beam direction of the particle beam by determining that the amount of the high-frequency component is equal to or less than the component amount threshold. The beam direction determination unit can determine that the direction of the beam source corresponding to the beam direction image of the electron density at which the zeroth-order moment is equal to or less than the moment threshold is an appropriate direction as the beam direction of the particle beam. Further, the beam direction determination unit can determine that the direction of the beam source corresponding to the beam direction image of the electron density at which the zero-order moment is larger than the moment threshold is an inappropriate direction as the particle beam direction.

  When the zeroth-order moment is used as an evaluation index, an appropriate beam direction can be determined with higher accuracy by simple calculation than when the absolute value of the average power spectrum slope is used as the evaluation index.

  The method for determining the beam direction of the particle beam using the computer in the treatment planning apparatus for creating the treatment plan information used in the particle beam irradiation system of the present invention implements the following first to third steps. In the first step, a predetermined angular interval is provided around the three-dimensional electron density image created based on the CT value obtained from the treatment plan CT image of the irradiation target obtained by imaging the irradiation target with the CT imaging apparatus. Assuming installation positions of a plurality of beam sources, the planned target is projected from the surface of the irradiation target in the three-dimensional electron density image on a two-dimensional screen virtual in the direction from the installation position of the beam source to the image of the planned target volume. A two-dimensional image obtained by projecting the electron density up to the end point of the volume is created for each installation position of a plurality of beam sources as a beam direction image of the electron density viewed from the installation position of the beam source. . In the second step, a power spectrum is calculated for each of the beam direction images having a plurality of electron densities. In the third step, the beam direction of the particle beam is determined based on the amount of the high frequency component in the spatial frequency component of the power spectrum. In the third step, the direction of the beam source capable of obtaining a beam direction image in which the amount of the high-frequency component is equal to or smaller than a predetermined component amount threshold is determined as an appropriate direction as the beam direction of the particle beam.

  When the absolute value of the average power spectrum slope is used as the evaluation index, in the second step, the average power spectrum is calculated for each of the beam direction images having a plurality of electron densities. Then, in the third step, the absolute value of the average power spectrum slope obtained for each of the plurality of electron density beam direction images is calculated, and the absolute value of the average power spectrum slope is greater than or equal to a predetermined slope threshold value. Determines that the amount of the high-frequency component is equal to or less than the component amount threshold, and determines the beam direction of the particle beam.

  When the 0th-order moment of the power spectrum is used as the evaluation index, in the second step, the integrated power spectrum is calculated for each of the beam direction images having a plurality of electron densities. In the third step, the zero-order moment is calculated from the integrated power spectrum obtained for each of the beam direction images having a plurality of electron densities, and when the zero-order moment is equal to or less than a predetermined zero-order moment threshold, Is determined to be equal to or less than the component amount threshold value, and the beam direction of the particle beam is determined.

  A computer program according to the present invention causes a computer to realize the first to third steps. This program can be stored in a storage medium.

(A) And (B) is a block diagram which shows the structure of an example of embodiment of the particle beam beam direction determination system of this invention used with the treatment plan apparatus which produces the treatment plan information used for a particle beam irradiation system. . FIG. 2 is a schematic diagram conceptually showing a situation in which an electron density beam direction image is created from a three-dimensional electron density image in order to facilitate understanding of the embodiment of FIG. 1. It is a flowchart which shows the algorithm of the software used when a beam direction image creation part acquires the beam direction image of an electron density. It is a figure which shows that the image of a plan target volume turns into a smooth isotropic image by CUBIC interpolation. It is a figure used in order to demonstrate obtaining the three-dimensional electron density image by sampling the electron density on the beam line from the surface of irradiation object to each terminal point of a plan target volume. The specific steps of the fourth step ST14 are shown in two dimensions. It is a figure used in order to explain the method of producing the beam direction image of electron density. (A) is a figure which shows the radiograph containing the image of the plan target volume of the head of the person who is irradiation object, (B) is a figure which shows the beam direction image of the electron density containing the image of the plan target volume. . It is a figure which shows the result obtained in an average power spectrum calculation part, an inclination calculation part, and a beam direction determination part in a time series, respectively. It is a flowchart which shows the algorithm of the software used when implement | achieving an average power spectrum calculation part and an inclination calculation part using a computer. (A) is the figure which substituted the average value of the electron density value in the image of the beam direction image of an electron beam to the image peripheral part of the beam direction image, (B) is a graph figure which shows the gray value in an image. is there. (A) And (B) is a figure which shows the change of the image and gray value before smoothing and after smoothing. (A) is a figure which shows the normalized image, (B) is a figure which shows the change of the gray value after the normalization before normalization. It is a figure which shows the replacement pattern of the image area | region of an original power spectrum image. It is a figure which shows a power spectrum image. It is a figure used in order to explain a power spectrum image. It is a figure used in order to explain a power spectrum image. (A) is a power spectrum image before conversion, and (B) is a diagram showing a result of polar coordinate conversion of the image of FIG. 18 (A). It is a figure which shows the result of having calculated | required the integral power spectrum from the power spectrum image which carried out the polar coordinate conversion shown in FIG.18 (B). It is a figure which shows the average power spectrum calculated | required from the integrated power spectrum. It is a figure which shows the algorithm of the software used by the method which determines the beam direction of a particle beam using a computer in the treatment plan apparatus which produces the treatment plan information used for a particle beam irradiation system. It is a figure which shows the gradient of the power spectrum with respect to the beam direction of 0-355 degree | times obtained in the test. (A) thru | or (C) are the figures which displayed the case where a beam was irradiated with the inclination of A thru | or C exceeding threshold value SH on FIG. 22 on three CT slices. It is a figure which shows an example of a structure of an irradiation field formation apparatus in broad beam irradiation. It is a figure which shows an example of the correspondence of the compensation filter and the plan target volume PTV which is a target. It is a block diagram which shows the structure of embodiment which determines whether the quantity of a high frequency component becomes below a predetermined component quantity threshold value based on a 0th moment. It is a figure which shows the procedure of the whole evaluation by a 0th moment. It is a figure which shows the detail of pre-processing. It is a figure which shows the algorithm of the conversion from a Cartesian coordinate system to a polar coordinate system. It is a figure which shows the power spectrum of the Cartesian coordinate system before conversion, and the power spectrum of the polar coordinate system after conversion. It is a figure which shows the power spectrum of the electron density beam direction image in 1 direction of 1 case. It is a figure which shows the power spectrum of the electron density beam direction image in two directions of 1 case. It is a figure which shows the 0th-order moment of the power spectrum from all directions in the 1st case. Comparison of the average value in the beam direction selected by the radiation oncologist and the avoided beam direction of the zero-order moment normalized to 0 and the maximum value in 5 cases to which the proposed method was applied It is a table | surface which shows. (A) And (B) is a figure which shows the relationship between the beam direction in the 1st case, and a primary moment and a secondary moment. It is a figure which shows the example of treatment plan CT imaged by tilting the body.

  Hereinafter, an example of an embodiment of a particle beam direction determination system and a particle beam direction determination method of the present invention will be described in detail with reference to the drawings. FIG. 1A is a block diagram showing a basic configuration of a particle beam beam direction determination system of the present invention used in a treatment plan apparatus for creating treatment plan information used in a particle beam irradiation system. In the basic configuration of FIG. 1A, the particle beam direction determination system includes a beam direction image creation unit 1, a power spectrum calculation unit 2, and a beam direction determination unit 3. The output of the beam direction determining means 3 is provided to the treatment planning device 4. The beam direction image creation unit 1 is created based on the CT value obtained from the treatment plan CT image of the irradiation target obtained by the three-dimensional electron density image generation system 5 imaging the irradiation target with the CT imaging device. Use density images. The power spectrum calculation unit 2 calculates a power spectrum for each of the plurality of beam direction images of the electron density viewed from the positions of the plurality of beam sources. In the present invention, an integrated power spectrum and / or an average power spectrum is used as the power spectrum. The present invention is that the power spectrum can be used to determine whether or not there is an electron density region portion that changes sharply with a large amplitude change in the beam direction image (whether or not the beam direction is robust against setup errors). It is based. Therefore, for example, when the arrangement position of the beam source is changed at an interval of 5 ° around the three-dimensional electron density image (if the incident direction of the beam is changed at an interval of 5 ° around the three-dimensional electron density image) ), An average power spectrum is calculated for each of 72 electron density beam direction images. The specific power spectrum calculation unit 2 first substitutes the average value of the electron density values in the image of the beam direction image of the electron density into the image peripheral part of the beam direction image, and smoothes the edge of the beam direction image. In the smoothing, a Gaussian filter is applied around the contour of the region having the value of the beam direction image to smooth the edge of the region having the value of the beam direction image. By this smoothing, it is prevented that an originally unnecessary high-frequency component is included in the average power spectrum from an edge generated when an electron density beam direction image is created. Next, preprocessing for normalizing the average value of the entire smoothed image to 0 is performed. Specifically, the average value of the entire image is subtracted from the value of the entire image. Then, a power spectrum image is obtained by calculation based on an image in which the average value is normalized to 0 (preprocessing is performed). The power spectrum image is imaged by performing a two-dimensional Fourier transform on the preprocessed image to obtain a power spectrum by calculation. In other words, the preprocessed electron density beam direction image data is Fourier transformed to obtain power for each frequency bandwidth (frequency resolution) and image it to obtain a power spectrum image. Then, the image areas of the power spectrum image are exchanged so that the zero frequency comes to the center of the image. Next, polar coordinate conversion of the power spectrum is performed to convert the power spectrum image into an image represented by coordinates of spatial frequency and angle. Next, an integrated power spectrum at each spatial frequency of the power spectrum image subjected to polar coordinate conversion is calculated. An average power spectrum is obtained from this integrated power spectrum. The average power spectrum is obtained by taking the natural logarithm of the vertical axis value and the horizontal axis value of the figure showing the integrated power spectrum. The beam direction determining means 3 determines the beam direction of the particle beam based on the amount of the high frequency component in the spatial frequency component of the power spectrum. That is, the beam direction determining means 3 determines that the direction of the beam source capable of obtaining a beam direction image in which the amount of the high frequency component is equal to or less than a predetermined component amount threshold is an appropriate direction as the beam direction of the particle beam. . Therefore, according to the present invention, it is possible to determine a beam direction that is robust against a setup error (a direction in which the dose distribution is less likely to collapse due to the setup error). Whether or not the amount of the high frequency component is equal to or less than a predetermined component amount threshold can be determined based on several evaluation indexes. One evaluation index is the absolute value of the slope of the average power spectrum, and the other evaluation index is a zeroth-order moment obtained from the integrated power spectrum.

  FIG. 1B is a block diagram showing a configuration of an embodiment for determining whether or not the amount of the high frequency component is equal to or less than a predetermined component amount threshold based on the absolute value of the slope of the average power spectrum. is there. FIG. 2 is a schematic diagram conceptually showing a situation where a beam direction image BDI of electron density is created from a three-dimensional electron density image in order to facilitate understanding of the embodiment of FIG. In FIG. 2, the three-dimensional shape of the irradiation object G is shown instead of the three-dimensional electron density image 3DI. Further, in FIG. 2, a target (tumor) TG is shown as corresponding to the planned target volume PTV. In the embodiment of FIG. 1, the particle beam direction determination system includes a beam direction image creation unit 1, a power spectrum calculation unit 2 including an average power spectrum calculation unit 2 ′, an inclination calculation unit 3A, and a beam direction determination unit 3B. The beam direction determining means 3 is composed of The beam direction image creation unit 1 is created based on the CT value obtained from the treatment plan CT image of the irradiation target obtained by the three-dimensional electron density image generation system 5 imaging the irradiation target with the CT imaging device. A density image 3DI is used.

  In FIG. 3, in the three-dimensional electron density image creation system 5, a software is used to create a three-dimensional electron density image 3DI using a computer, and the beam direction image creation unit 1 acquires a beam direction image BDI of electron density. 2 is a flowchart showing a wear algorithm. Conceptually, the three-dimensional electron density image creation system 5 is a three-dimensional electron density image based on a CT value obtained from a treatment plan CT image of an irradiation target obtained by imaging the irradiation target with a CT imaging device (not shown). You get 3DI. In the present embodiment, a three-dimensional electron density image 3DI is specifically created by the following steps. As shown in FIG. 3, in the first step ST11, isotropic voxelization of the treatment plan CT image and the plan target volume in the treatment plan CT image is performed. In particular, at the time of isotropic voxelization, as shown in FIG. 4, the image of the planned target volume PTV that has been isotropically voxeled is converted into a smooth isotropic image by CUBIC interpolation. This suppresses an increase in high-frequency components. Next, in a second step, the CT value of the treatment plan CT image converted into isotropic voxels is converted into the relative electron density of water to create a converted electron density image. The conversion to electron density is specifically disclosed, for example, on pages 250 and 251 of “Radiation Therapy Physics” (author: Takehiro Nishijo) issued by Bunkodo Co., Ltd. In the third step, a three-dimensional electron density image 3DI is obtained by extracting a body region (non-air region) portion from the three-dimensional conversion electron density image. In step ST4, as shown in FIG. 5, the range from the surface of the irradiation target to each terminal point of the planned target volume PTV is calculated using the three-dimensional electron density image 3DI, and based on this range. Then, the electron density on the beam line obtained by sampling is added to create a beam direction image of the electron density.

  FIG. 6 shows the specific steps of the fourth step ST14 in two dimensions. In FIG. 6, B is an image of the bed, and W is an image of a net for fixing the patient to the bed. These images B and W are erased by image processing. Specifically, fixed threshold processing of converted electron density image → opening process (1 pixel reduction / 1 pixel enlargement) → maximum area cutout → closing process (1 pixel enlargement / 1 pixel reduction) → image inversion and maximum area cutout (outside of body) Air) → “region including air in the body” is defined as a body region → the body region is extracted from the electron density image. When the three-dimensional electron density image 3DI is created through such processing, a three-dimensional electron density image 3DI that more accurately reflects the electron density of the beam passing portion can be obtained.

  As shown conceptually in FIG. 7, the beam direction image creation unit 1 in FIG. 1B provides a predetermined angular interval (5 ° interval in the present embodiment) around the three-dimensional electron density image 3DI. Assuming the installation positions (P0 to P71) of the plurality of beam sources LS, the three-dimensional electron density on the two-dimensional screen S virtually assumed in the direction from the installation position of each beam source LS toward the image PTVI of the planned target volume PTV. A two-dimensional image obtained by projecting the electron density from the surface of the irradiation target in the image 3DI to the end point of the planned target volume PTV is a beam direction image BDI of the electron density viewed from the position of the beam source. Create for each installation position of multiple beam sources. In FIGS. 7 and 2, F corresponds to the multi-leaf collimator shown in FIG. 24, and is for projecting a beam corresponding to the contour shape of the planned target volume PTV onto the three-dimensional electron density image 3DI. The two-dimensional screen S is a virtual calculation screen for obtaining a beam direction image by calculation, and the projection of the three-dimensional electron density image 3DI onto the two-dimensional screen is also a calculation projection. In the beam direction image BDI of the electron density, the electron density in the three-dimensional electron density image on the beam line from the beam source LS side is represented by the terminal point of the planned target volume PTV from the surface of the irradiation target (the planned target volume viewed from the beam source). A value obtained by adding the samples sampled up to the point corresponding to the rear edge of the PTV) is included as two-dimensional data for each of a large number of beam lines. If there is a cavity in the region through which the beam passes, the sum of the electron density sampled on the beam passing through the region including the cavity in the electron direction beam direction image BDI will be low, and regions such as muscle and bone will be removed. The sum of the electron densities sampled on the passing beam is high. Incidentally, FIG. 8A shows a radiograph including an image of the planned target volume of the head of the person to be irradiated, and FIG. 8B shows an electron density beam direction image BDI including the image of the planned target volume. It is. The beam direction image creation unit 1 acquires such a beam direction image BDI having an electron density by changing the installation position of the beam source LS from 0 ° to 355 ° at intervals of 5 °.

  The average power spectrum calculation unit 2 ′ in FIG. 1 calculates an average power spectrum for each of the beam direction images BDI having a plurality of electron densities viewed from the arrangement positions of the plurality of beam sources LS in FIG. As shown in FIG. 7, assuming the position of the beam source at 5 ° intervals around the three-dimensional electron density image 3DI, the average power spectrum is calculated for each of 72 electron density beam direction images. Become. The inclination calculation unit 3A calculates the inclination of the average power spectrum obtained from the beam direction images having a plurality of electron densities. Then, the beam direction determination unit 3B determines the beam direction of the particle beam based on the absolute value of this inclination. The information determined by the beam direction determination unit 3B is provided to the treatment planning device 4.

  FIG. 9 shows the results obtained by the average power spectrum calculation unit 2 ′, the inclination calculation unit 3A, and the beam direction determination unit 3B in time series. FIG. 10 is a flowchart showing a software algorithm used when the average power spectrum calculation unit 2 ′ and the inclination calculation unit 3A are realized using a computer. The upper diagram in FIG. 11A shows a diagram in which the average value of the electron density values in the image of the beam direction image of the electron beam is substituted into the image peripheral portion of the beam direction image. In addition, the horizontal line in a figure shows the part which extract | collected the gray value mentioned later. Next, in step ST21, the average power spectrum calculation unit 2 'substitutes the average value of the electron density values in the image of the beam direction image of the beam into the image peripheral part of the beam direction image BDI, and the edge of the beam direction image BDI is substituted. Smooth. This reduces the difference between the pixel values of the background image around the beam direction image BDI and the pixel values of the beam direction image BDI, so that the beam direction image based on the electron density inside the planned target volume PTV. This means that the average pixel value of BDI is assigned to the pixel value of the background image. The result of step ST21 can be confirmed as a change in gray value in the image, as shown in FIG. In FIG. 11B, the gray value in the image changes from the state of the upper graph to the state of the lower graph.

  Next, in step ST22, the edge of the projection part is smoothed. In this smoothing, a Gaussian filter is applied around the contour of the region having the value of the beam direction image BDI to smooth the edge of the region having the value of the beam direction image BDI. That is, a large value difference is eliminated between the portion having the value of the beam direction image BDI and the portion not having the value, thereby preventing an increase in the high frequency component. Here, the Gaussian filter is used, for example, to reduce high frequency components in the power spectrum in a predetermined pixel area from the end of the beam direction image BDI of electron density. FIGS. 12A and 12B show changes in the image and gray value before and after smoothing.

  Next, in step ST23, preprocessing for normalizing the average value of the entire smoothed image to 0 is performed. Specifically, the average value of the entire image is subtracted from the value of the entire image. FIG. 13A shows a normalized image, and FIG. 13B shows a change in gray value after normalization before normalization.

  Next, in step ST24, a power spectrum image is obtained by calculation on the basis of an image in which the average value is normalized to 0 (preprocessing is performed). The power spectrum image is obtained by performing a two-dimensional Fourier transform on the preprocessed image to obtain a power spectrum by calculation and imaging it. In other words, the preprocessed electron density beam direction image data is Fourier transformed to obtain power for each frequency bandwidth (frequency resolution) and image it to obtain a power spectrum image. In the present embodiment, the image areas are switched so that the power spectrum image has 0 frequency at the center of the image. FIG. 14 shows a replacement pattern of image areas of the original power spectrum image. In the original power spectrum image, the four corners have zero frequency. Therefore, the region R1 of the image is replaced with the region R4, and the region R2 and the region R3 are replaced. As a result, the zero frequency is located at the center of the image. FIG. 15 shows an example of a power spectrum image obtained by exchanging image areas. The density fluctuation of the power spectrum image obtained in this way is indicated by changes in frequency and amplitude. The coordinates (distance from the center) of the power spectrum image indicate the spatial frequency (frequency of density fluctuation of the image). A small power spectrum means that the amplitude is small. As shown in FIG. 15, low frequency components gather at the center of the power spectrum image, and high frequency components gather at the outer region. In FIG. 16, the spatial frequency of the region indicated by the symbol A has a low frequency component, and the spatial frequency of the region indicated by the symbol B is a high frequency component. Further, as in the region indicated by reference numeral C in FIG. 17, a portion having a large value at the same distance from the center indicates that the density change is large. In addition, as in the region indicated by the symbol D, a portion having a small value at the same distance from the center indicates that the density change is small.

Next, in step ST25, polar coordinate conversion of the power spectrum image is performed to convert the power spectrum image into an image represented by coordinates of spatial frequency and angle. FIG. 18A shows a power spectrum image before conversion, and FIG. 18B shows the result of polar coordinate conversion of the image of FIG. In step ST26, an integrated power spectrum at each spatial frequency of the power spectrum image subjected to polar coordinate conversion is calculated. FIG. 19 shows the result of obtaining the integrated power spectrum from the polar spectrum converted power spectrum image shown in FIG. The integrated power spectrum can be obtained by the following equation.

  In the above formula, f is a spatial frequency, and θ is an angle.

  Step ST27 is a step for realizing the inclination calculation unit 3A. In step ST27, the slope of the average power spectrum is calculated from the integrated power spectrum obtained from the beam direction images of a plurality of electron densities. The average power spectrum is obtained by taking the natural logarithm of the vertical axis value and the horizontal axis value of the figure showing the integrated power spectrum. The slope or slope of the average power spectrum is obtained as the slope of a first order polynomial with the power spectrum for all frequencies up to the Nyquist frequency. FIG. 20 shows an average power spectrum obtained from the integrated power spectrum of FIG. The absolute value of the slope of the average power spectrum is obtained as the slope of the approximate line SL by obtaining an approximate line SL from the average power spectrum by the least square method.

  The beam direction determination unit 3B determines the beam direction of the particle beam based on the absolute value of this inclination. The lower the spatial frequency of the change in the beam direction image of electron density and the smaller the amplitude, the larger the absolute value of the gradient or slope of the average power spectrum. When the region of the electron density that changes sharply with a large amplitude change in the beam direction image becomes smaller, the spatial frequency of the change in the electron direction beam direction image becomes lower and the amplitude becomes smaller. The high frequency component in the spatial frequency component is reduced, and the absolute value of the gradient, that is, the slope of the average power spectrum is increased. Therefore, if the absolute value of this slope is large, even if there is some patient setup error, there is no region of electron density that changes sharply with large amplitude changes in the beam direction image. The likelihood of occurrence is reduced.

  Therefore, the beam direction determination unit 3B determines that the direction of the beam source corresponding to the beam direction image of the electron density at which the absolute value of the tilt is larger than the tilt threshold SH is an appropriate direction as the beam direction of the particle beam. According to the method of determining the tilt threshold, a plurality of types of beam directions are determined as preferable directions. Further, the beam direction determination unit 3B may determine that the direction of the beam source corresponding to the beam direction image having an electron density smaller than the inclination threshold SH is an inappropriate direction as the particle beam direction. it can. If it does in this way, what is necessary is just to employ | adopt the direction which a doctor considers preferable from directions other than an inappropriate direction.

  FIG. 21 shows a software algorithm used in a method for determining the beam direction of a particle beam using a computer in a treatment plan apparatus for creating treatment plan information used in the particle beam irradiation system of the present invention. In this algorithm, in a first step ST31, a predetermined area around a three-dimensional electron density image created based on a CT value obtained from a treatment plan CT image of an irradiation target obtained by imaging the irradiation target with a CT imaging apparatus. Assuming the installation positions of a plurality of beam sources with an angular interval of, the irradiation target in the three-dimensional electron density image is placed on a virtual two-dimensional screen in the direction from the installation position of each beam source to the planned target volume. A two-dimensional image obtained by projecting the electron density from the surface to the end point of the planned target volume is created for each of a plurality of beam sources as a beam direction image of the electron density viewed from the position of the beam source. . In the second step ST32, an average power spectrum is calculated for each of the beam direction images having a plurality of electron densities. In the third step ST33, the beam direction of the particle beam is determined based on the absolute value of the gradient of the average power spectrum obtained from the beam direction images having a plurality of electron densities.

  In order to confirm the effect of the present invention, the present embodiment was applied to five head and neck cancer cases, and a test for determining the beam direction was performed. FIG. 22 shows the power spectrum slope for the 0-355 degree beam direction obtained in this test. FIGS. 23A to 23C are diagrams showing the case where the beams are irradiated with the inclinations A to C exceeding the threshold SH in FIG. 22 on three CT slices. When we asked the radiation oncologist for their views on these three beam directions, they agreed that it was preferable. This means that the beam direction determined by the present invention indicates a robust beam direction.

  FIG. 26 is a block diagram showing a configuration of an embodiment for determining whether or not the amount of the high-frequency component is equal to or less than a predetermined component amount threshold based on the 0th-order moment. When the 0th-order moment is used as the evaluation index, the power spectrum calculation unit 2 calculates a one-dimensional integrated power spectrum for each of a plurality of density beam direction images. The beam direction determining means 30 includes a 0th moment calculation unit 30A for calculating a 0th moment from the one-dimensional integral power spectrum obtained for each of the beam direction images having a plurality of electron densities, and the 0th moment is below a predetermined moment threshold. In this case, it is determined that the amount of the high-frequency component is equal to or less than the component amount threshold, and the beam direction determining unit 30B determines the beam direction of the particle beam. The beam direction determination unit 30B can determine that the direction of the beam source corresponding to the beam direction image of the electron density at which the zero-order moment is equal to or less than the moment threshold is an appropriate direction as the beam direction of the particle beam. Further, the beam direction determination unit 30B can determine that the direction of the beam source corresponding to the beam direction image of the electron density at which the zero-order moment is greater than the moment threshold is an inappropriate direction as the beam direction of the particle beam. .

  Hereinafter, a method for determining the beam direction using the 0th-order moment will be described in detail. As described above, the degree of dose distribution deterioration due to patient setup error increases with the amplitude and spatial frequency of the spatial variation of the water equivalent distance. Therefore, in the present embodiment, a beam direction having an electron density beam direction image having a low spatial frequency component and a small amplitude is found using the zeroth-order moment. The zero-order moment (power spectrum area) of the one-dimensional integrated power spectrum (non-patent document 16) obtained from the electron density beam direction image in each beam direction is obtained. FIG. 27 shows the overall procedure of evaluation using the zeroth moment. First, preprocessing is performed on the acquired electron density beam direction image. Next, the two-dimensional power spectrum of the Cartesian coordinate system is calculated. Thereafter, the power spectrum of the Cartesian coordinate system is converted to polar coordinates. Then, the integrated power spectrum is calculated, and finally the zeroth moment is calculated.

  First, the power spectrum of the electron density beam direction image was calculated from 0 to 355 degrees at intervals of 5 degrees. Before calculating the power spectrum, the following pre-processing is performed on the electron density beam direction image in order to reduce high-frequency components that are unrelated to fluctuations in the water equivalent distance in the boundary region of the virtual irradiation field. Details of the preprocessing are shown in FIG. First, in order to reduce the difference between the values in the virtual irradiation field and the outside of the irradiation field, the average pixel value in the virtual irradiation field is substituted into the virtual irradiation field. Next, a Gaussian filter was applied to 10 pixels inside and outside the virtual irradiation field. In order to reduce the influence of the average value of the image on the power spectrum, normalization was performed by subtracting the average pixel value of the entire image from the entire image.

A two-dimensional power spectrum was calculated by applying a two-dimensional Fourier transform to the preprocessed image. The two-dimensional Fourier transform and power spectrum were calculated according to the following formula.

  Here, f (x, y) is a pre-processed electron density beam direction image, F (u, v) is a Fourier transform image, P (u, v) is a power spectrum image, and x and y are in real space. The coordinates of u and v are coordinates in the frequency space.

Next, the power spectrum P (f, θ) of the polar coordinate system was obtained using the following equation.

  In general, since the coordinates of (u, v) are real numbers, the value of P (u, v) was obtained using a linear interpolation method. FIG. 29 shows an algorithm for conversion from a Cartesian coordinate system to a polar coordinate system, and FIG. 30 shows a power spectrum of a Cartesian coordinate system before conversion and a power spectrum of a polar coordinate system after conversion.

From the two-dimensional power spectrum P (f, θ), the angle was integrated from 0 to 360 degrees, and the integrated power spectrum P (f) was calculated. The calculation formula is as follows.

Here, f is the spatial frequency (mm ⁻¹ ) in the polar coordinate power spectrum, and θ is the direction (°). Further, since the original power spectrum has a very large value near zero frequency and the value of the high frequency component is relatively very small, the natural logarithm of the original power spectrum was used.

Finally, the zero-order moment M ₀ was calculated from the one-dimensional integrated power spectrum (Non-patent Document 16). The calculation of the zero-order moment M ₀ is described in detail in “Gerzberg L, Meindl JD 1980 Power-spectrum centroid detection for doppler systems applications Ultrasonic Imaging 2 232-261”. The calculation formula is as follows.

  In order to confirm the effect of this embodiment, this embodiment was applied to five cases of nasal cavity tumor. First, the power spectrum of the electron density beam direction image was compared between the beam direction chosen by the radiation oncologist and the avoidable beam direction. In general, radiation oncologists plan their treatment by avoiding beams that pass through the nasal septum in parallel or through the cones of the temporal bone. 31 and 32 show the power spectrum of the electron density beam direction image in the two directions in the first case. In the figure, a CT image, a PTV, a path in the beam direction selected by the radiation oncologist, and a path in the beam direction to be avoided are shown. In FIG. 31, the power spectrum in the 0 degree beam direction parallel to the nasal septum, which is avoided by the radiation oncologist, especially in the high frequency region, is larger than the 35 degree direction value chosen by the radiation oncologist. . This result means that the electron density beam direction image in the 0 degree direction includes many high-frequency components, and irradiation from the 0 degree direction means that the dose distribution is likely to be deformed by a patient setup error. Also, in FIG. 32, the radiation oncologist selected the 115 degree direction rather than the beam direction passing through the cone like the 160 degree direction. As in the case of FIG. 31, the power spectrum in the 160-degree beam direction avoided by the radiation oncologist is larger than the value in the selected beam direction in the high-frequency region. The results shown in FIGS. 31 and 32 show that evaluation using the zero-order moment of the power spectrum is useful for quantitative evaluation of the beam direction.

  FIG. 33 shows the zeroth-order moment of the power spectrum from all directions in the first case. As shown in the figure, the value of the beam direction chosen by the radiation oncologist (35, 115, 250, 335 degrees direction indicated by the dotted line) is small compared to other beam directions that are avoided by the radiation oncologist. It is a value.

  In the table of FIG. 34, in the 5 cases to which the proposed method is applied, the zero-order moment normalized to 0 and the maximum value of 1 are normalized with the beam direction selected by the radiation oncologist and the avoidable beam direction. A comparison of mean values is shown.

  The other two indexes are included in the table of FIG. 34, which will be described later. Other indices are normalized in the same manner as the zeroth-order moment. As a result, the zero-order moment is statistically the average value in the beam direction selected by the radiation oncologist between the average value in the beam direction selected by the radiation oncologist and the average value in the beam direction avoided. It was significantly smaller. Therefore, it can be said that this embodiment can quantitatively select a beam direction that is robust against patient setup errors based on the relationship between the beam direction and the zero-order moment of the power spectrum of the electron density beam direction image.

  As shown in the table of FIG. 34, there was no statistically significant difference in the average value of the first and second moments between the selected beam direction and the avoided beam direction. FIGS. 35A and 35B show the relationship between the beam direction, the first moment, and the second moment in the first case. Unlike the 0th moment, there was no characteristic tendency in the 1st and 2nd moments, and the values oscillated. Therefore, it can be said that only the 0th-order moment is suitable for selecting a beam direction that is robust to patient setup errors.

  In carbon beam therapy, there are many cases where there is no rotating gantry, and the beam irradiation ports are only horizontal, vertical, and 45 degrees. Therefore, before the treatment plan, the beam direction is determined empirically by the planner and the body is tilted. After creating the fixture, the treatment plan CT must be imaged. An example of a treatment plan CT imaged while tilting the body is shown in FIG. After imaging the treatment plan CT, if it is determined that the beam direction is inappropriate, the fixture must be created again and the treatment plan CT must be re-imaged. However, patients undergoing particle beam therapy usually have a diagnostic CT examination in advance. By applying this embodiment to the diagnostic CT, the beam direction can be quantitatively determined before imaging the treatment plan CT.

  In the above embodiment, the 0th order moment is obtained from the one-dimensional integrated power spectrum, but it is needless to say that the 0th order moment of the power spectrum may be calculated from the two-dimensional integrated power spectrum.

  According to the present invention, it is possible to provide a particle beam beam direction determination system, a particle beam beam direction determination method, and a particle beam beam direction determination computer program that determine a beam direction that is robust against patient setup errors. .

DESCRIPTION OF SYMBOLS 1 Beam direction image creation part 2 Power spectrum calculation part 2 'Average power spectrum calculation part 3 Beam direction determination means 3A Inclination calculation part 3B Beam direction determination part 4 Treatment planning apparatus 5 Three-dimensional electron density image creation system

Claims (19)

A particle beam beam direction determination system for determining a beam direction of a particle beam, which is used in a treatment planning apparatus for creating treatment plan information used for a particle beam irradiation system,
A plurality of beam sources at predetermined angular intervals around a three-dimensional electron density image created based on a CT value obtained from a treatment plan CT image of the irradiation target obtained by imaging the irradiation target with a CT imaging apparatus On the two-dimensional screen virtual in the direction from the installation position of the beam source to the image of the planned target volume, the planned target volume from the surface of the irradiation target in the three-dimensional electron density image A two-dimensional image obtained by projecting the electron density up to the end point of the beam is created for each installation position of the plurality of beam sources as a beam direction image of the electron density viewed from the installation position of the beam source. A beam direction image creation unit,
A power spectrum calculation unit for calculating a power spectrum for each of the plurality of beam direction images of the electron density,
A particle beam beam direction determination system comprising beam direction determination means for determining a beam direction of the particle beam based on an amount of a high frequency component in a spatial frequency component of the power spectrum.   The beam direction determining means is such that the direction of the beam source capable of obtaining a beam direction image in which the amount of the high-frequency component is equal to or less than a predetermined component amount threshold is an appropriate direction as the beam direction of the particle beam. The particle beam beam direction determination system according to claim 1 for determining. The power spectrum calculation unit is composed of an average power spectrum calculation unit that calculates an average power spectrum for each of the plurality of beam direction images of the electron density,
The beam direction determining means includes
A slope calculating unit for calculating an absolute value of a slope of an average power spectrum obtained for each of a plurality of beam direction images of the electron density;
Beam direction for determining the beam direction of the particle beam by determining that the amount of the high-frequency component is equal to or less than the component amount threshold when the absolute value of the gradient of the average power spectrum is equal to or greater than a predetermined gradient threshold The particle beam beam direction determination system according to claim 2, comprising a determination unit. The beam direction determination unit includes:
The direction of the beam source corresponding to the beam direction image of the electron density at which the absolute value of the tilt is equal to or greater than the tilt threshold is determined as an appropriate direction as the beam direction of the particle beam. 4. The particle beam direction determining system according to 3. The beam direction determination unit includes:
The direction of the beam source corresponding to the beam direction image of the electron density whose absolute value of the tilt is smaller than the tilt threshold is determined to be an inappropriate direction as the beam direction of the particle beam. Item 4. The particle beam direction determination system according to Item 3. The power spectrum calculation unit calculates an integrated power spectrum for each of the plurality of electron density beam direction images,
The beam direction determining means includes
A 0th-order moment calculation unit for calculating a 0th-order moment from the integrated power spectrum obtained for each of the plurality of beam densities of the electron density;
When the zero-order moment is equal to or less than a predetermined moment threshold, the beam direction determination unit determines that the amount of the high-frequency component is equal to or less than the component amount threshold and determines the beam direction of the particle beam. Item 3. The particle beam direction determination system according to Item 2. The beam direction determination unit includes:
The direction of the beam source corresponding to the beam direction image of the electron density at which the zero-order moment is equal to or less than the moment threshold is determined as an appropriate direction as the beam direction of the particle beam. A particle beam beam direction determination system as described in. The beam direction determination unit includes:
The direction of the beam source corresponding to the beam direction image of the electron density at which the zero-order moment is greater than the moment threshold is determined to be an inappropriate direction as the beam direction of the particle beam. Item 7. The beam direction determining system according to Item 6. The three-dimensional electron density image is
Performing isotropic voxelization of the treatment plan CT image and the plan target volume in the treatment plan CT image;
Converting the CT value of the treatment plan CT image into the isotropic voxel into a relative electron density of water to create a converted electron density image,
The particle beam beam direction determination system according to any one of claims 1 to 8, wherein a body region is extracted from the converted electron density image.   The projection of the electron density in the beam direction image creation unit is executed by sampling and adding the electron density on a beam line from the surface of the irradiation target to each end point of the planned target volume. The particle beam beam direction determination system according to any one of 1 to 8. The average power spectrum calculation unit
Substituting an average value of electron density values in the beam direction image of the electron density into the image peripheral portion in the beam direction image, smoothing the edge of the beam direction image,
Normalize the average value of the entire smoothed image to 0,
Based on the image whose average is normalized to 0, a power spectrum image is obtained by calculation,
The coordinates of the power spectrum image are converted into polar coordinates,
Calculate the integrated power spectrum at each spatial frequency of the power spectrum image converted to polar coordinates,
The particle beam beam direction determination system according to claim 3, wherein the average power spectrum is obtained from the integrated power spectrum. A particle beam beam direction determining method for determining a beam direction of a particle beam, which is used in a treatment planning apparatus for generating treatment plan information used in a particle beam irradiation system,
A plurality of beam sources at predetermined angular intervals around a three-dimensional electron density image created based on a CT value obtained from a treatment plan CT image of the irradiation target obtained by imaging the irradiation target with a CT imaging apparatus On the two-dimensional screen virtual in the direction from the installation position of the beam source to the image of the planned target volume, the planned target volume from the surface of the irradiation target in the three-dimensional electron density image A two-dimensional image obtained by projecting the electron density up to the end point of the beam is created for each installation position of the plurality of beam sources as a beam direction image of the electron density viewed from the installation position of the beam source. A first step to:
A second step of calculating a power spectrum for each of a plurality of beam direction images of the electron density;
A particle beam direction determination method comprising: a third step of determining a beam direction of the particle beam based on an amount of a high frequency component in a spatial frequency component of the power spectrum.   In the third step, the direction of the beam source capable of obtaining a beam direction image in which the amount of the high-frequency component is equal to or less than a predetermined component amount threshold is an appropriate direction as the beam direction of the particle beam. 13. The particle beam direction determination method according to claim 12, wherein the particle beam direction is determined. In the second step, an average power spectrum is calculated for each of the plurality of electron density beam direction images,
In the third step,
Calculate the absolute value of the slope of the average power spectrum obtained for each of the plurality of electron density beam direction images,
The beam direction of the particle beam is determined by determining that the amount of the high-frequency component is equal to or less than the component amount threshold when the absolute value of the gradient of the average power spectrum is equal to or greater than a predetermined gradient threshold. 14. The particle beam direction determination method according to 13. In the second step, an integrated power spectrum is calculated for each of the plurality of beam direction images of the electron density,
In the third step, a zero-order moment is calculated from the integrated power spectrum obtained for each of a plurality of beam direction images of the electron density, and when the zero-order moment is equal to or less than a predetermined moment threshold, The particle beam direction determination method according to claim 13, wherein the beam direction of the particle beam is determined by determining that the amount of the high-frequency component is equal to or less than the component amount threshold value. The three-dimensional electron density image is
Performing isotropic voxelization of the treatment plan CT image and the plan target volume in the treatment plan CT image;
Converting the CT value of the treatment plan CT image into the isotropic voxel into a relative electron density of water to create a converted electron density image,
The particle beam direction determination method according to any one of claims 12 to 15, wherein a body region is extracted from the converted electron density image.   13. The projection of the electron density in the first step is performed by sampling and adding the electron density on a beam line from the surface of the irradiation object to each end point of the planned target volume. 16. The particle beam direction determination method according to any one of items 1 to 15. In the second step,
Substituting the average value of the electron density values in the image including the beam direction image of the electron density into the image peripheral portion of the beam direction image, smoothing the edge of the beam direction image,
Normalize the average value of the entire smoothed image to 0,
Based on the image whose average is normalized to 0, a power spectrum image is obtained by calculation,
The coordinates of the power spectrum image are converted into polar coordinates,
Calculate the integrated power spectrum at each spatial frequency of the power spectrum image converted to polar coordinates,
Obtain the average power spectrum from the integrated power spectrum,
The particle beam direction determination method according to claim 14, wherein an inclination of the average power spectrum is calculated. It is a program installed in a computer when a particle beam direction determination method for determining the beam direction of a particle beam, which is used in a treatment plan apparatus for creating treatment plan information used in a particle beam irradiation system, is executed using a computer. And
A plurality of beam sources at predetermined angular intervals around a three-dimensional electron density image created based on a CT value obtained from a treatment plan CT image of the irradiation target obtained by imaging the irradiation target with a CT imaging apparatus On the two-dimensional screen virtual in the direction from the installation position of the beam source to the image of the planned target volume, the planned target volume from the surface of the irradiation target in the three-dimensional electron density image A two-dimensional image obtained by projecting the electron density up to the end point of the beam is created for each installation position of the plurality of beam sources as a beam direction image of the electron density viewed from the installation position of the beam source. A first step to:
A second step of calculating a power spectrum for each of a plurality of beam direction images of the electron density;
A particle beam configured to cause a computer to execute a third step of determining a beam direction of the particle beam based on an amount of a high frequency component in a spatial frequency component of the power spectrum. Direction determining program.