CN109754860B - Method for automatically predicting planning difficulty level based on DVH prediction model - Google Patents
Method for automatically predicting planning difficulty level based on DVH prediction model Download PDFInfo
- Publication number
- CN109754860B CN109754860B CN201811596187.0A CN201811596187A CN109754860B CN 109754860 B CN109754860 B CN 109754860B CN 201811596187 A CN201811596187 A CN 201811596187A CN 109754860 B CN109754860 B CN 109754860B
- Authority
- CN
- China
- Prior art keywords
- risk
- pca
- curve
- organ
- dvh
- 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
Landscapes
- Measuring And Recording Apparatus For Diagnosis (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
The invention discloses a method for automatically predicting planning difficulty degree based on a DVH prediction model, which comprises the following steps: (1) collecting delineation data and planning data of the same type of cancer patients, and training a DVH prediction model of organs at risk by using the collected data; (2) inputting target areas of new cases and delineation data of organs at risk, and predicting DVH curves of the organs at risk by using a trained DVH prediction model; (3) the ease of planning is assessed in terms of the predicted organ-at-risk DVH curves. The method can effectively automate the process of predicting the difficulty level of the plan, thereby realizing automation of radiotherapy plan allocation and improving the accuracy and efficiency of predicting the difficulty level of the plan.
Description
Technical Field
The invention belongs to the field of radiotherapy, and particularly relates to a method for automatically predicting planning difficulty degree based on a DVH prediction model.
Background
The planning design of radiotherapy (hereinafter referred to as "radiotherapy") is an important part of the radiotherapy process, physical physicians with different experiences can design radiotherapy plans with different time consumption and different planning quality, and in actual work, in order to improve the efficiency and quality of the whole radiotherapy process, hospitals can distribute the plans with different difficulties to the designs of the physical physicians with different experiences, so that a problem naturally occurs: how to predict the ease of planning?
In the past radiotherapy process, the prediction of planning difficulty level is implemented by physicists based on own experience, and the traditional method has two disadvantages:
1. the workload of the physicist is increased;
2. the experience of the physicist is heavily dependent.
Disclosure of Invention
In order to solve the technical problems, the invention provides a method for automatically predicting the difficulty level of a plan based on a DVH prediction model, which can effectively automate the process of predicting the difficulty level of the plan, thereby realizing automation of radiotherapy plan allocation and improving the accuracy and efficiency of prediction of the difficulty level of the plan.
In order to achieve the purpose, the technical scheme of the invention is as follows:
a method for automatically predicting planning difficulty level based on DVH prediction model comprises the following steps:
(1) collecting delineation data and planning data of the same type of cancer patients, and training a DVH prediction model of organs at risk by using the collected data;
(2) inputting target areas of new cases and delineation data of organs at risk, and predicting DVH curves of the organs at risk by using a trained DVH prediction model;
(3) the ease of planning is assessed in terms of the predicted organ-at-risk DVH curves.
On the basis of the technical scheme, the following improvements can be made:
as a preferred scheme, in step (1), the training of the DVH prediction model specifically includes the following steps:
(1.1) calculating or acquiring OVH curves, DVH curves and volumes of all organs at risk of all patients, and acquiring the volumes of target areas of all patients;
(1.2) selecting OVH curves of all patients on the organs at risk for each organ at risk, training a PCA dimension reduction model of the OVH on the organs at risk by taking the OVH curves as training data, and simultaneously performing PCA dimension reduction on the OVH by using the trained model to obtain PCA characteristics of the OVH;
(1.3) training a PCA dimension reduction model of the DVH of each organ at risk in the same way as the step (1.2), and calculating PCA characteristics of the DVH of all patients on each organ at risk;
(1.4) for m patients, obtaining m training samples;
for each patient, a training sample was obtained as follows:
(1.4.1) sequentially arranging the volume of each organ at risk, the volume of a target area and each PCA characteristic of each organ at risk OVH of a patient to obtain one datum;
(1.4.2) sequentially arranging each PCA characteristic of each organ at risk DVH of the patient to obtain a label corresponding to the data;
(1.5) putting m training samples into the SVR model for training.
As a preferable scheme, the step (2) specifically comprises the following steps:
(2.1) inputting the sketching data of a patient;
(2.2) calculating OVH curves of each organ at risk of the patient, and performing PCA dimension reduction on each OVH curve by using a corresponding PCA dimension reduction model obtained in a DVH prediction model training stage to obtain PCA features of each OVH curve;
(2.3) acquiring the volume of each endangered organ and the target area of the patient;
wherein the calculation method of each volume is required to be consistent with the calculation method of the volume in the training stage of the DVH prediction model;
(2.4) sequentially arranging the volume of each organ at risk, the volume of a target area and each PCA characteristic of each organ at risk OVH of a patient to obtain data, and respectively inputting the data into each PCA characteristic prediction model of each organ at risk DVH to obtain a predicted value;
(2.5) for each organ at risk, sequentially arranging the predicted values of the PCA features of the organ at risk DVH to form a PCA feature vector of the organ at risk, and performing an inverse transformation of the PCA dimensionality reduction of the organ at risk DVH by taking the feature vector as an input to predict the DVH curve of the organ at risk.
Preferably, the method for calculating the patient organ-at-risk OVH curve comprises the following steps:
s1: interpolating the delineation into three-dimensional volume data:
s2: OVH curves were calculated on three-dimensional volume data.
Preferably, step S1 specifically includes the following steps:
s1.1: unifying the coordinate units of all the delineation points of the organs at risk and the target PTV into mm, namely multiplying the coordinate in the direction of the delineation point X, Y by the pixel of CT;
s1.2: traversing all sketches of the organs at risk and the target PTV, and respectively finding out minimum and maximum slice positions minP and maxP, minimum and maximum X coordinates minX and maxX, and minimum and maximum Y coordinates minY and maxY of the sketches;
s1.3: calculating the coordinate span range x-minX and range y-minY of the X, Y direction, and recording the maximum value of range x and range y as range;
s1.4: specifying a maximum pixel length pixelLen of the voxel in the direction X, Y, and then calculating therefrom a pixel size pixelSize of the voxel range/pixelLen;
s1.5: the pixel length of the three-dimensional volume in the X, Y, Z slice direction is calculated,
xLen＝ceil(rangeX/pixelSize)；
yLen＝ceil(rangeY/pixelSize)；
zLen＝ceil((maxP-minP)/pixelSize)；
s1.6: initializing a three-dimensional volume of the target region PTV to a zero matrix of scale (xLen, yLen, zLen);
s1.7: traversing all sketches in the target region PTV, calculating a Z coordinate ptz corresponding to the slice position where the sketches are located in the three-dimensional body, simultaneously calculating pixel points of all sketching points in the sketches in an XOY plane corresponding to the three-dimensional body and marking the pixel points as pts, and filling polygons in a two-dimensional image of which the Z coordinate of the three-dimensional body of the target region PTV is the ptz by using pts, wherein the filling value is 1;
if (x, y) is a delineation point on the delineation with the slice position p, the calculation formula of the corresponding pixel point (ptx, pty, ptz) is:
ptx＝round((x-minX)/pixelSize)；
pty＝round((y-minY)/pixelSize)；
ptz＝round((p-minP)/ppixelSize)；
s1.8: interpolating a three-dimensional volume of the target PTV in the Z direction, specifically:
acquiring Z coordinates minZ and maxZ corresponding to the minimum and maximum slice positions delineated by the target region PTV in a three-dimensional body, interpolating a two-dimensional graph of which the Z coordinates between the minZ and the maxZ are all 0 in the three-dimensional body, and selecting the value of a two-dimensional image corresponding to a slice closest to the two-dimensional image in the Z direction as the value of the two-dimensional image;
s1.9: steps S1.6-S1.8 are also performed on the organs at risk, and their corresponding three-dimensional volumes are obtained.
Preferably, step S2 specifically includes the following steps:
s2.1: acquiring the boundary of the PTV in the target area three-dimensional body by adopting a six-neighborhood method, specifically, when a pixel point in the PTV in the target area is on the boundary of the three-dimensional body or points with the value of 0 exist in the upper part, the lower part, the front part, the rear part and the left part and the right part of the PTV, the PTV is considered as a boundary point;
s2.2: calculating the distance from each point of the three-dimensional body of the organ at risk to the PTV boundary of the target area, and specifically comprising the following steps:
s2.2.1: obtaining an internal part inPart and an external part outPart of a endangered organ in a three-dimensional body in a target region PTV;
s2.2.2: respectively calculating the distance from each point in InPart and outPart to the boundary of the target region PTV, wherein the distance of each point in InPart is expressed in the PTV by taking the inverse number;
s2.3: arranging the distances in an ascending order, acquiring minimum maximum values minD and maxD in the distances, and then acquiring t equal division points d0, d1, d2, a.
S2.4: finding the subscript ki of the first element larger than di in the ascending distance list, wherein the subscript starts from 0, ki is 100/pixelNum is the ordinate corresponding to di in the OVH curve, and the subscript ki is the subscript ki of the first element larger than di, wherein: pixelNum is the number of pixels of the organ at risk in the three-dimensional volume; when i is taken over 0, 1, 2.. times.t, an OVH curve of the organ at risk can be obtained.
As a preferred scheme, in the step (1.2) and the step (1.3), a PCA dimension reduction model is trained, and the trained PCA dimension reduction model is used to respectively perform dimension reduction on an OVH curve and a DVH curve to obtain corresponding PCA characteristics, specifically comprising the following steps:
t1: extracting initial characteristics of the curve, uniformly adopting a fixed number of sample points on the curve, and discretizing the curve;
t2: training a curve PCA dimension reduction model;
t3: and (4) executing PCA dimension reduction of the curve, extracting initial characteristics of the curve, and then carrying out PCA dimension reduction on the initial characteristics by using parameters in the PCA dimension reduction model.
Preferably, step T1 specifically includes the following steps:
t1.1: only the point with the maximum corresponding abscissa is reserved in the points with repeated ordinate, and other points are removed from the curve;
t1.2: the ordinate is normalized to [0, A ] by the following method, wherein minY and maxY are respectively the minimum value and the maximum value of the ordinate,
y：＝(y-minY)/(maxY-minY)；
t1.3: and uniformly taking a numbers in [0, A ] to form a vertical coordinate vector (0, 1A/(a-1), 2A/(a-1).., A), and linearly interpolating to obtain a corresponding horizontal coordinate vector of the vertical coordinate vector in the curve, wherein the horizontal coordinate vector is called as the initial characteristic of the curve.
Preferably, step T2 specifically includes the following steps:
t2.1: if m curves participate in training, respectively extracting initial features of the m curves to obtain m initial feature vectors which are respectively marked as X1, X2,. and Xm;
t2.2: calculating a mean vector X _ mean and a standard deviation vector X _ std of the m initial feature vectors, wherein the ith component of the X _ std is the standard deviation of a vector formed by the ith components of the m initial feature vectors;
t2.3: carrying out standardization processing on the initial feature vector;
Xi：＝(Xi-X_mean)/X_std，i＝1，2，..，m；
wherein, the division of the above formula refers to the division of the corresponding elements;
t2.4: calculating a sample covariance matrix C (cij) n of the m eigenvectors;
wherein, cij is cov (Xi, Xj), i, j is 1, 2.
T2.5: calculating eigenvalues and eigenvectors of the covariance matrix C, selecting the eigenvalues with the largest absolute values, wherein the ratio of the sum of the absolute values of the eigenvalues to the sum of the absolute values of all eigenvalues is greater than a threshold;
the unit eigenvectors corresponding to these eigenvalues are called principal components of the PCA, and are respectively denoted as p1, p 2.
Preferably, step T3 specifically includes the following steps:
t3.1: inputting a curve, extracting the initial characteristics of the curve and marking as X;
t3.2: the following operations are performed for X as follows,
X：＝(X-X_mean)/X_std；
wherein, the division of the above formula is the division of the corresponding components, and X _ mean and X _ std are obtained in the model training process of the step T2;
t3.3: if X is a row vector, PCA principal components p1, p2 and k obtained in the model training process are column vectors, the projection of X on each principal component is calculated,
X：＝X*(p1，p2，...，pk)＝(X*p1，X*p2，...，X*pk)；
wherein: x is the characteristic vector after the dimensionality reduction of curve PCA and is called the PCA characteristic vector of the curve, each component of X is called the PCA characteristic of the curve, and the multiplication in the above formula refers to matrix multiplication.
Drawings
Fig. 1 is a flowchart of a method for automatically predicting a planning difficulty level based on a DVH prediction model according to an embodiment of the present invention.
Fig. 2 is a flowchart of DVH prediction model training according to an embodiment of the present invention.
Fig. 3 is a flowchart of a DVH curve prediction method for an organ at risk according to an embodiment of the present invention.
Fig. 4 is a flow chart illustrating the ease of planning a DVH assessment based on organ-at-risk prediction according to an embodiment of the present invention.
Fig. 5 is a scoring criteria chart for prostate planning provided by an embodiment of the present invention.
Fig. 6 is a flowchart of curve initial feature extraction according to an embodiment of the present invention.
Fig. 7 is a flowchart of training a curve PCA dimension reduction model according to an embodiment of the present invention.
Fig. 8 is a flow chart of the dimensionality reduction of the curve PCA according to the embodiment of the present invention.
Fig. 9 is a flowchart of the inverse dimension reduction of the curve PCA according to the embodiment of the present invention.
Detailed Description
Preferred embodiments of the present invention will be described in detail below with reference to the accompanying drawings.
To achieve the object of the present invention, in some embodiments of a method for automatically predicting difficulty of planning based on a DVH prediction model, as shown in fig. 1, the method for automatically predicting difficulty of planning based on a DVH prediction model includes the following steps:
(1) collecting delineation data and planning data of the same type of cancer patients, and training a DVH prediction model of organs at risk by using the collected data;
(2) inputting target areas of new cases and delineation data of organs at risk, and predicting DVH curves of the organs at risk by using a trained DVH prediction model;
(3) the ease of planning is assessed in terms of the predicted organ-at-risk DVH curves.
In step (1), to collect intensity modulated plans for a sufficient number of patients with the same type of cancer (which ensures that the patients have the same organs at risk), one plan is collected for each case, and to prevent errors in field angles, the field angles of the plans are required to be the same, and to ensure the quality of the training data, each plan needs to be designed by or approved by experienced physicists.
In step (1), a Support Vector Regression (SVR) is used to calculate the relationship between PCA features of the organ at risk DVH and patient anatomical features, where the anatomical features include: volume of each organ at risk, target volume, PCA signature of each organ at risk OVH. As shown in fig. 2, the training of the DVH prediction model specifically includes the following steps:
(1.1) calculating or acquiring OVH curves, DVH curves and volumes of all organs at risk of all patients, and acquiring the volumes of target areas of all patients;
(1.2) selecting OVH curves of all patients on the organs at risk for each organ at risk, training a PCA dimension reduction model of the OVH on the organs at risk by taking the OVH curves as training data, and simultaneously performing PCA dimension reduction on the OVH by using the trained model to obtain PCA characteristics of the OVH;
(1.3) training a PCA dimension reduction model of the DVH of each organ at risk in the same way as the step (1.2), and calculating PCA characteristics of the DVH of all patients on each organ at risk;
(1.4) for m patients, obtaining m training samples;
for each patient, a training sample was obtained as follows:
(1.4.1) sequentially arranging the volume of each organ at risk, the volume of a target area and each PCA characteristic of each organ at risk OVH of a patient to obtain one datum;
(1.4.2) sequentially arranging each PCA characteristic of each organ at risk DVH of the patient to obtain a label corresponding to the data;
(1.5) putting m training samples into the SVR model for training.
And (3) calling an existing Support Vector Regression (SVR) model in sklern of python by the SVR model related to the step (1.5), and in addition, because the SVR in sklern requires that the output is a scalar, splitting the DVH prediction model into a plurality of submodels for training, wherein the number of the submodels is the sum of the PCA characteristic numbers of all organs at risk DVH, the input of the submodels is unchanged, and the output is the PCA characteristics of each organ at risk DVH.
So far, three models can be obtained in the DVH prediction model training process: a PCA dimension reduction model of the OVH, a PCA dimension reduction model of the DVH, and a PCA feature prediction model of the DVH, which together form a DVH prediction model.
In step (2), the execution of the DVH prediction model calculates various anatomical features of the patient according to the patient delineation, and then calculates a DVH curve of the organs at risk of the patient using the trained SVR relationship model, as shown in fig. 3, which specifically includes the following steps:
(2.1) inputting the sketching data of a patient;
(2.2) calculating OVH curves of each organ at risk of the patient, and performing PCA dimension reduction on each OVH curve by using a corresponding PCA dimension reduction model obtained in a DVH prediction model training stage to obtain PCA features of each OVH curve;
(2.3) acquiring the volume of each endangered organ and the target area of the patient;
wherein the calculation method of each volume is required to be consistent with the calculation method of the volume in the training stage of the DVH prediction model;
(2.4) sequentially arranging the volume of each organ at risk, the volume of a target area and each PCA characteristic of each organ at risk OVH of a patient to obtain data, and respectively inputting the data into each PCA characteristic prediction model of each organ at risk DVH to obtain a predicted value;
(2.5) for each organ at risk, sequentially arranging the predicted values of the PCA features of the organ at risk DVH to form a PCA feature vector of the organ at risk, and performing an inverse transformation of the PCA dimensionality reduction of the organ at risk DVH by taking the feature vector as an input to predict the DVH curve of the organ at risk.
In step (3), the difficulty level of the plan is evaluated according to the predicted DVH of the organs at risk, a scoring standard is made according to the plan competition method, the predicted DVH of the organs at risk is scored according to the scoring standard, the higher the score is, the easier the plan is considered to be designed, and conversely, the plan is considered to be difficult to be designed. As shown in fig. 4, the specific steps are as follows:
(3.1) making a scoring criterion for the plan following the plan tournament, fig. 5 being the scoring criterion for the prostate plan we made;
and (3.2) calculating scores of all DVH indexes of all organs at risk according to scoring standards, and then summing the scores to obtain the total score of all current organs at risk DVH.
Here, taking the DVH index V75 of Rectum (Rectum) in the prostate plan as an example, a method for calculating the score is described,
if V75 is less than 10%, the score is 5;
if V75 > 15%, score 0;
if V75 is in the interval [ 10%, 15% ], using a linear interpolation method to obtain the score;
(3.3) normalizing score to 100 according to the maximum possible value of score, comprising the following steps:
(3.3.1) calculate score max _ score maximum possible (max _ score value 5 x 8-40 in prostate plan if only two organs at risk, Rectum (recovery) and Bladder (Bladder), are considered);
(3.3.2)score：＝score*100/max_score；
(3.4) divide the ease of planning into five grades by size of score: very difficult, moderate, easy, very easy, and thus assess the ease of planning;
the division criteria are as follows:
it is very difficult: score has a value of 0;
it is difficult to: score is greater than 0 and less than or equal to 33.3;
the method is moderate: score is greater than 33.3 and less than or equal to 66.6;
the method is easy to realize: score size 66.6, less than 100;
it is very easy: score has a value of 100.
For a better understanding of the present invention, the present invention is further described below, including but not limited to the following:
the definition of patient organs at risk OVH (Overlap Volume Histogram) is as follows:
in the formula:
t is a target area;
o is an organ at risk;
| O | is the volume of the organ at risk;
p is a subset of O;
d (p, T) is the distance of p to the tumor boundary;
{ p ∈ O | d (p, T) ≦ T } representing a set of voxels in the organ at risk O that are less than distance T from target T;
the overlapping volume histogram function of the target T and the organ at risk O is the volume fraction of the distance from the organ at risk O to the target T less than T.
The method for calculating the OVH curve of the organs at risk of the patient comprises the following steps:
s1: interpolating the delineation into three-dimensional volume data:
s2: OVH curves were calculated on three-dimensional volume data.
Step S1 specifically includes the following steps:
s1.1: unifying the coordinate units of all the delineation points of the organs at risk and the target PTV into mm, namely multiplying the coordinate in the direction of the delineation point X, Y by the pixel of CT;
s1.2: traversing all sketches of the organs at risk and the target PTV, and respectively finding out minimum and maximum slice positions minP and maxP, minimum and maximum X coordinates minX and maxX, and minimum and maximum Y coordinates minY and maxY of the sketches;
s1.3: calculating the coordinate span range x-minX and range y-minY of the X, Y direction, and recording the maximum value of range x and range y as range;
s1.4: specifying the maximum pixel length pixelLen of the voxel in the direction X, Y (generally taking 128), and then calculating therefrom the pixel size pixelSize range/pixelLen of the voxel;
s1.5: the pixel length of the three-dimensional volume at X, Y, Z (slice direction) is calculated,
xLen＝ceil(rangeX/pixelSize)；
yLen＝ceil(rangeY/pixelSize)；
zLen＝ceil((maxP-minP)/pixelSize)；
the function ceil returns the smallest integer greater than or equal to the specified expression;
s1.6: initializing a three-dimensional volume of the target region PTV to a zero matrix of scale (xLen, yLen, zLen);
s1.7: traversing all sketches in the target region PTV, calculating a Z coordinate ptz corresponding to the slice position where the sketches are located in the three-dimensional body, simultaneously calculating pixel points of all sketching points in the sketches in an XOY plane corresponding to the three-dimensional body and marking the pixel points as pts, and filling polygons in a two-dimensional image of which the Z coordinate of the three-dimensional body of the target region PTV is the ptz by using pts, wherein the filling value is 1;
if (x, y) is a delineation point on the delineation with the slice position p, the calculation formula of the corresponding pixel point (ptx, pty, ptz) is:
ptx＝round((x-minX)/pixelSize)；
pty＝round((y-minY)/pixelSize)；
ptz＝round((p-minP)/ppixelSize)；
the function round functions to return a rounded value for a specified expression;
s1.8: interpolating a three-dimensional volume of the target PTV in the Z direction, specifically:
acquiring Z coordinates minZ and maxZ corresponding to the minimum and maximum slice positions delineated by the target region PTV in a three-dimensional body, interpolating a two-dimensional graph of which the Z coordinates between the minZ and the maxZ are all 0 in the three-dimensional body, and selecting the value of a two-dimensional image corresponding to a slice closest to the two-dimensional image in the Z direction as the value of the two-dimensional image;
s1.9: steps S1.6-S1.8 are also performed on organs at risk (e.g., rectal rectum, bladder blader) to obtain their corresponding three-dimensional volumes.
Step S2 specifically includes the following steps:
s2.1: acquiring the boundary of the PTV in the target area three-dimensional body by adopting a six-neighborhood method, specifically, when a pixel point in the PTV in the target area is on the boundary of the three-dimensional body or points with the value of 0 exist in the upper part, the lower part, the front part, the rear part and the left part and the right part of the PTV, the PTV is considered as a boundary point;
s2.2: calculating the distance from each point of the three-dimensional body of the organ at risk to the PTV boundary of the target area, and specifically comprising the following steps:
s2.2.1: obtaining an internal part inPart and an external part outPart of a endangered organ in a three-dimensional body in a target region PTV;
s2.2.2: respectively calculating the distance from each point in InPart and outPart to the boundary of the target region PTV, wherein the distance of each point in InPart is expressed in the PTV by taking the inverse number;
s2.3: arranging the distances in an ascending order, acquiring minimum maximum values minD and maxD in the distances, and then acquiring 50 equally dividing points d0, d1, d2, and d50 of an interval [ minD, maxD ];
s2.4: finding the subscript ki (subscript starts from 0) of the first element larger than di in the ascending distance list, wherein ki x 100/pixelNum is the ordinate corresponding to di in the OVH curve, and wherein: pixelNum is the number of pixels of the organ at risk in the three-dimensional volume; when i is taken over 0, 1, 2.., 50, an OVH curve of the organ at risk is obtained.
In the step (1.2) and the step (1.3), a PCA dimension reduction model is trained, and the trained PCA dimension reduction model is used for respectively reducing dimensions of an OVH curve and a DVH curve to obtain corresponding PCA characteristics, and the method specifically comprises the following steps:
t1: extracting initial characteristics of the curve, uniformly adopting a fixed number of sample points on the curve, and discretizing the curve;
t2: training a curve PCA dimension reduction model;
t3: and (4) executing PCA dimension reduction of the curve, extracting initial characteristics of the curve, and then carrying out PCA dimension reduction on the initial characteristics by using parameters in the PCA dimension reduction model.
As shown in fig. 6, step T1 specifically includes the following:
t1.1: only the point with the maximum corresponding abscissa is reserved in the points with repeated ordinate, and other points are removed from the curve;
t1.2: the ordinate is normalized to [0, 100] by the following method, wherein minY and maxY are respectively the minimum value and the maximum value of the ordinate,
y：＝(y-minY)/(maxY-minY)；
t1.3: taking 50 numbers uniformly in [0, 100] to form a vertical coordinate vector (0, 1 × 100/49, 2 × 100/49., 100), and linearly interpolating to obtain a corresponding horizontal coordinate vector of the vertical coordinate vector in the curve, wherein the horizontal coordinate vector is called as an initial feature of the curve, namely a 50-dimensional vector.
As shown in fig. 7, step T2 specifically includes the following:
t2.1: if m curves participate in training, respectively extracting initial features of the m curves to obtain m initial feature vectors which are respectively marked as X1, X2,. and Xm;
t2.2: calculating a mean vector X _ mean and a standard deviation vector X _ std of the m initial feature vectors, wherein the ith component of the X _ std is the standard deviation of a vector formed by the ith components of the m initial feature vectors;
t2.3: carrying out standardization processing on the initial feature vector;
Xi：＝(Xi-X_mean)/X_std，i＝1，2，..，m；
wherein, the division of the above formula refers to the division of the corresponding elements;
t2.4: calculating a sample covariance matrix C (cij) n of the m eigenvectors;
wherein, cij is cov (Xi, Xj), i, j is 1, 2.
T2.5: calculating eigenvalues and eigenvectors of the covariance matrix C, selecting the eigenvalues with the largest absolute values, wherein the ratio of the sum of the absolute values of the eigenvalues to the sum of the absolute values of all eigenvalues is greater than a threshold;
the unit eigenvectors corresponding to these eigenvalues are called principal components of the PCA, and are respectively denoted as p1, p 2.
So far, the training of the curve PCA dimension reduction model is completed, and the mean vector X _ mean and the standard deviation vector X _ std of the initial feature vector calculated in step T2.2 and the PCA principal components p1, p 2.
The threshold value threshold in step T2.5 is specified externally, and its size determines the number k of PCA principal components, typically 0.95. In addition, the value of the number k of principal components of PCA can be directly specified during training.
As shown in fig. 8, step T3 specifically includes the following:
t3.1: inputting a curve, extracting the initial characteristics of the curve and marking as X;
t3.2: the following operations are performed for X as follows,
X：＝(X-X_mean)/X_std；
wherein, the division of the above formula is the division of the corresponding components, and X _ mean and X _ std are obtained in the model training process of the step T2;
t3.3: if X is a row vector, PCA principal components p1, p2 and so, which are obtained in the model training process, and pk is a column vector (the assumed conditions are not met and corresponding adjustment can be made), calculating the projection of X on each principal component,
X：＝X*(p1，p2，...，pk)＝(X*p1，X*p2，...，X*pk)；
wherein: x is the characteristic vector after the dimensionality reduction of curve PCA and is called the PCA characteristic vector of the curve, each component of X is called the PCA characteristic of the curve, and the multiplication in the above formula refers to matrix multiplication.
The inverse transformation of the dimensionality reduction of the curve PCA is described below, and the overall idea is to restore the PCA characteristics to initial characteristics by using parameters of a PCA dimensionality reduction model, then obtain a curve by using the initial characteristics as abscissa and using sampling points of the ordinate when the curve initial characteristics are acquired as ordinate, as shown in fig. 9, the specific steps are as follows:
f1: a PCA feature vector X is input (pcl, pc 2.., pck), X is restored to the initial feature space,
X：＝pc1*p1+pc2*p2+...+pck*pk；
X：＝(XT*X_std)+X_mean；
the multiplication of the above formula is corresponding component multiplication, XT is transfer operation, X _ mean and X _ std are mean vector and standard deviation vector obtained by PCA dimension reduction model respectively (assuming row vector, not adjusted correspondingly), p1, p 2.
F2: taking X as an abscissa vector, (0, 1 × 100/49, 2 × 100/49.,. 100) as an ordinate vector (the vector is the ordinate vector when the initial features of the curve are collected), a curve can be obtained, and this is the restored curve.
In summary, the present invention constructs an SVR relationship model between patient anatomical features and patient organ-at-risk DVH (dose volume histogram) and uses this model to predict the patient organ-at-risk DVH, assigns scoring criteria for the plan in the manner of a plan competition, and scores the predicted DVH, with higher scores considering the plan as easier to design and vice versa considering the plan as harder to design. The invention can realize the process automation of predicting the difficulty degree of the plan, thereby realizing the automation of the distribution of the radiotherapy plan and improving the accuracy and efficiency of predicting the difficulty degree of the plan.
With respect to the preferred embodiments of the present invention, it should be noted that, for those skilled in the art, various changes and modifications can be made without departing from the inventive concept of the present invention, and these changes and modifications are within the scope of the present invention.
Claims (5)
1. A method for automatically predicting planning difficulty level based on DVH prediction model is characterized by comprising the following steps:
(1) collecting delineation data and planning data of the same type of cancer patients, and training a DVH prediction model of organs at risk by using the collected data;
(2) inputting target areas of new cases and delineation data of organs at risk, and predicting DVH curves of the organs at risk by using a trained DVH prediction model;
(3) evaluating the difficulty of planning according to the predicted organ-at-risk DVH curve;
in the step (1), the training of the DVH prediction model specifically includes the following steps:
(1.1) calculating or acquiring OVH curves, DVH curves and volumes of all organs at risk of all patients, and acquiring the volumes of target areas of all patients;
(1.2) selecting OVH curves of all patients on the organs at risk for each organ at risk, training a PCA dimension reduction model of the OVH on the organs at risk by taking the OVH curves as training data, and simultaneously performing PCA dimension reduction on the OVH by using the trained model to obtain PCA characteristics of the OVH;
(1.3) training a PCA dimension reduction model of the DVH of each organ at risk in the same way as the step (1.2), and calculating PCA characteristics of the DVH of all patients on each organ at risk;
(1.4) for m patients, obtaining m training samples;
for each patient, a training sample was obtained as follows:
(1.4.1) sequentially arranging the volume of each organ at risk, the volume of a target area and each PCA characteristic of each organ at risk OVH of a patient to obtain one datum;
(1.4.2) sequentially arranging each PCA characteristic of each organ at risk DVH of the patient to obtain a label corresponding to the data;
(1.5) putting the m training samples into an SVR model for training;
the step (2) specifically comprises the following steps:
(2.1) inputting the sketching data of a patient;
(2.2) calculating OVH curves of each organ at risk of the patient, and performing PCA dimension reduction on each OVH curve by using a corresponding PCA dimension reduction model obtained in a DVH prediction model training stage to obtain PCA features of each OVH curve;
(2.3) acquiring the volume of each endangered organ and the target area of the patient;
wherein the calculation method of each volume is required to be consistent with the calculation method of the volume in the training stage of the DVH prediction model;
(2.4) sequentially arranging the volume of each organ at risk, the volume of a target area and each PCA characteristic of each organ at risk OVH of a patient to obtain data, and respectively inputting the data into each PCA characteristic prediction model of each organ at risk DVH to obtain a predicted value;
(2.5) for each organ at risk, sequentially arranging the predicted values of the PCA features of the organ at risk DVH to form a PCA feature vector of the organ at risk, and performing inverse transformation of PCA dimensionality reduction of the organ at risk DVH by taking the feature vector as input to predict a DVH curve of the organ at risk;
the method for calculating the OVH curve of the organs at risk of the patient comprises the following steps:
s1: interpolating the delineation into three-dimensional volume data:
s2: calculating an OVH curve on the three-dimensional volume data;
the step S1 specifically includes the following steps:
s1.1: unifying the coordinate units of all the delineation points of the organs at risk and the target PTV into mm, namely multiplying the coordinate in the direction of the delineation point X, Y by the pixel of CT;
s1.2: traversing all sketches of the organs at risk and the target PTV, and respectively finding out minimum and maximum slice positions minP and maxP, minimum and maximum X coordinates minX and maxX, and minimum and maximum Y coordinates minY and maxY of the sketches;
s1.3: calculating the coordinate span range X = maxX-minX and range Y = maxY-minY of the direction X, Y, and recording the maximum value of range X and range Y as range;
s1.4: specifying a maximum pixel length pixelLen of the voxel in the direction X, Y, and then calculating therefrom a pixel size pixelSize = range/pixelLen of the voxel;
s1.5: the pixel length of the three-dimensional volume in the X, Y, Z slice direction is calculated,
xLen=ceil(rangeX/pixelSize)；
yLen=ceil(rangeY/pixelSize)；
zLen=ceil(( maxP-minP)/ pixelSize)；
s1.6: initializing a three-dimensional volume of the target region PTV to a zero matrix of scale (xLen, yLen, zLen);
s1.7: traversing all sketches in the target region PTV, calculating a Z coordinate ptz corresponding to the slice position where the sketches are located in the three-dimensional body, simultaneously calculating pixel points of all sketching points in the sketches in an XOY plane corresponding to the three-dimensional body and marking the pixel points as pts, and filling polygons in a two-dimensional image of which the Z coordinate of the three-dimensional body of the target region PTV is the ptz by using pts, wherein the filling value is 1;
if (x, y) is a delineation point on the delineation with the slice position p, the calculation formula of the corresponding pixel point (ptx, pty, ptz) is:
ptx=round((x-minX)/ pixelSize)；
pty=round((y-minY)/ pixelSize)；
ptz=round((p-minP)/ppixelSize)；
s1.8: interpolating a three-dimensional volume of the target PTV in the Z direction, specifically:
acquiring Z coordinates minZ and maxZ corresponding to the minimum and maximum slice positions delineated by the target region PTV in a three-dimensional body, interpolating a two-dimensional graph of which the Z coordinates between the minZ and the maxZ are all 0 in the three-dimensional body, and selecting the value of a two-dimensional image corresponding to a slice closest to the two-dimensional image in the Z direction as the value of the two-dimensional image;
s1.9: the same steps S1.6-S1.8 are carried out on the organs at risk, and corresponding three-dimensional bodies are obtained;
the step S2 specifically includes the following steps:
s2.1: acquiring the boundary of the PTV in the target area three-dimensional body by adopting a six-neighborhood method, specifically, when a pixel point in the PTV in the target area is on the boundary of the three-dimensional body or points with the value of 0 exist in the upper part, the lower part, the front part, the rear part and the left part and the right part of the PTV, the PTV is considered as a boundary point;
s2.2: calculating the distance from each point of the three-dimensional body of the organ at risk to the PTV boundary of the target area, and specifically comprising the following steps:
s2.2.1: obtaining an internal part inPart and an external part outPart of a endangered organ in a three-dimensional body in a target region PTV;
s2.2.2: respectively calculating the distance from each point in InPart and outPart to the boundary of the target region PTV, wherein the distance of each point in InPart is expressed in the PTV by taking the inverse number;
s2.3: arranging the distances in an ascending order, obtaining the minimum maximum values minD and maxD in the distances, and then obtaining t equal division points d0, d1, d2, … and dt of an interval [ minD, maxD ];
s2.4: finding the subscript ki of the first element larger than di in the ascending distance list, wherein the subscript starts from 0, ki is 100/pixelNum is the ordinate corresponding to di in the OVH curve, and the subscript ki is the subscript ki of the first element larger than di, wherein: pixelNum is the number of pixels of the organ at risk in the three-dimensional volume; when i is taken 0, 1, 2, …, t, an OVH curve of the organ at risk is obtained.
2. The method for automatically predicting planning difficulty based on the DVH prediction model according to claim 1, wherein in the step (1.2) and the step (1.3), a PCA dimension reduction model is trained, and an OVH curve and a DVH curve are respectively dimension reduced by using the trained PCA dimension reduction model to obtain corresponding PCA features, specifically comprising the following steps:
t1, extracting initial characteristics of the curve, uniformly taking a fixed number of sample points on the curve, and discretizing the curve;
t2, training a curve PCA dimension reduction model;
and T3, executing PCA dimension reduction of the curve, extracting initial characteristics of the curve, and then carrying out PCA dimension reduction on the initial characteristics by using parameters in the PCA dimension reduction model.
3. The method according to claim 2, wherein said step T1 comprises the following steps:
t1.1: only the point with the maximum corresponding abscissa is reserved in the points with repeated ordinate, and other points are removed from the curve;
t1.2: the ordinate is normalized to [0, A ] by the following method, wherein minY and maxY are respectively the minimum value and the maximum value of the ordinate,
y: = (y-minY)/(maxY-minY)；
t1.3: and uniformly taking a number in [0, A ] to form a vertical coordinate vector (0, 1A/(a-1), 2A/(a-1) …, A), and linearly interpolating to obtain a corresponding horizontal coordinate vector of the vertical coordinate vector in the curve, wherein the horizontal coordinate vector is called as the initial characteristic of the curve.
4. The method according to claim 2, wherein said step T2 comprises the following steps:
t2.1: if m curves participate in training, respectively extracting initial features of the m curves to obtain m initial feature vectors which are respectively marked as X1, X2, … and Xm;
t2.2: calculating a mean vector X _ mean and a standard deviation vector X _ std of the m initial feature vectors, wherein the ith component of the X _ std is the standard deviation of a vector formed by the ith components of the m initial feature vectors;
t2.3: carrying out standardization processing on the initial feature vector;
Xi := (Xi - X_mean)/X_std，i=1, 2, .., m；
wherein, the division of the above formula refers to the division of the corresponding elements;
t2.4: calculating a sample covariance matrix C = (cij) n = n of m eigenvectors;
wherein cij = cov (Xi, Xj), i, j =1, 2, …, n;
t2.5: calculating eigenvalues and eigenvectors of the covariance matrix C, selecting the eigenvalues with the largest absolute values, wherein the ratio of the sum of the absolute values of the eigenvalues to the sum of the absolute values of all eigenvalues is greater than a threshold;
the unit eigenvectors corresponding to these eigenvalues are called principal components of the PCA, and are respectively denoted as p1, p 2.
5. The method according to claim 2, wherein said step T3 comprises the following steps:
t3.1, inputting a curve, extracting the initial characteristics of the curve and marking as X;
t3.2 for X,
X := (X - X_mean)/ X_std；
wherein, the division of the above formula is the division of the corresponding components, and X _ mean and X _ std are obtained in the model training process of the step T2;
t3.3, if X is a row vector, PCA principal components p1, p2, … and pk obtained in the model training process are column vectors, calculating the projection of X on each principal component,
X := X*(p1, p2, …, pk) = (X*p1, X*p2, …, X*pk)；
wherein: x is the characteristic vector after the dimensionality reduction of curve PCA and is called the PCA characteristic vector of the curve, each component of X is called the PCA characteristic of the curve, and the multiplication in the above formula refers to matrix multiplication.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811596187.0A CN109754860B (en) | 2018-12-21 | 2018-12-21 | Method for automatically predicting planning difficulty level based on DVH prediction model |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811596187.0A CN109754860B (en) | 2018-12-21 | 2018-12-21 | Method for automatically predicting planning difficulty level based on DVH prediction model |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109754860A CN109754860A (en) | 2019-05-14 |
| CN109754860B true CN109754860B (en) | 2020-11-10 |
Family
ID=66404120
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811596187.0A Active CN109754860B (en) | 2018-12-21 | 2018-12-21 | Method for automatically predicting planning difficulty level based on DVH prediction model |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109754860B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110232964B (en) * | 2019-06-04 | 2023-11-14 | 苏州雷泰智能科技有限公司 | Plan implementation method and device based on predicted dose guidance and Gaussian process optimization |
| CN110503043A (en) * | 2019-08-24 | 2019-11-26 | 东莞市强艺体育器材有限公司 | Noise reduction homing method and system based on depth camera acquisition human body movement data |
| CN114904164A (en) * | 2022-05-10 | 2022-08-16 | 苏州雷泰医疗科技有限公司 | Radiotherapy automatic plan generation method and device based on artificial intelligence |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105031833A (en) * | 2015-08-28 | 2015-11-11 | 瑞地玛医学科技有限公司 | Dosage verification system for radiotherapy apparatus |
| CN105358219A (en) * | 2013-05-21 | 2016-02-24 | 瓦里安医疗系统国际股份公司 | Systemes and methods for automatic creation of dose prediction models and therapy treatment plans as a cloud service |
| CN107403201A (en) * | 2017-08-11 | 2017-11-28 | 强深智能医疗科技（昆山）有限公司 | Tumour radiotherapy target area and jeopardize that organ is intelligent, automation delineation method |
| CN107441637A (en) * | 2017-08-30 | 2017-12-08 | 南方医科大学 | The intensity modulated radiation therapy Forecasting Methodology of 3-dimensional dose distribution and its application in the works |
| CN108629785A (en) * | 2018-05-10 | 2018-10-09 | 西安电子科技大学 | Based on the three-dimensional magnetic resonance pancreas image partition method from step study |
| CN108766563A (en) * | 2018-05-25 | 2018-11-06 | 戴建荣 | Radiotherapy prediction of result method and system based on dosage group |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8688618B2 (en) * | 2009-06-23 | 2014-04-01 | The Johns Hopkins University | Method and system for determining treatment plans |
| EP3395405B1 (en) * | 2017-04-24 | 2021-10-13 | RaySearch Laboratories AB | System and method for radiation therapy treatment planning |
-
2018
- 2018-12-21 CN CN201811596187.0A patent/CN109754860B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105358219A (en) * | 2013-05-21 | 2016-02-24 | 瓦里安医疗系统国际股份公司 | Systemes and methods for automatic creation of dose prediction models and therapy treatment plans as a cloud service |
| CN105031833A (en) * | 2015-08-28 | 2015-11-11 | 瑞地玛医学科技有限公司 | Dosage verification system for radiotherapy apparatus |
| CN107403201A (en) * | 2017-08-11 | 2017-11-28 | 强深智能医疗科技（昆山）有限公司 | Tumour radiotherapy target area and jeopardize that organ is intelligent, automation delineation method |
| CN107441637A (en) * | 2017-08-30 | 2017-12-08 | 南方医科大学 | The intensity modulated radiation therapy Forecasting Methodology of 3-dimensional dose distribution and its application in the works |
| CN108629785A (en) * | 2018-05-10 | 2018-10-09 | 西安电子科技大学 | Based on the three-dimensional magnetic resonance pancreas image partition method from step study |
| CN108766563A (en) * | 2018-05-25 | 2018-11-06 | 戴建荣 | Radiotherapy prediction of result method and system based on dosage group |
Non-Patent Citations (3)
| Title |
|---|
| "A planning quality evaluation tool for prostate adaptive IMRT based on machine learning";ZHU X, GE Y, LI T, et al;《Medical Physics》;20111230;第719-726页 * |
| "基于神经网络学习方法的放疗计划三维剂量分布预测";孔繁图等;《南方医科大学学报》;20180627;第683-690页 * |
| 孔繁图等."基于神经网络学习方法的放疗计划三维剂量分布预测".《南方医科大学学报》.2018, * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109754860A (en) | 2019-05-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Cardenas et al. | Advances in auto-segmentation | |
| JP6514325B2 (en) | System and method for segmenting medical images based on anatomical landmark-based features | |
| EP3961565A1 (en) | Evaluating quality of segmentation of an image into different types of tissue for planning treatment using tumor treating fields (ttfields) | |
| CN109754860B (en) | Method for automatically predicting planning difficulty level based on DVH prediction model | |
| CN106920234B (en) | Combined automatic radiotherapy planning method | |
| EP3052188B1 (en) | Predicting achievable dose distribution using 3d information as an input | |
| CN104424629B (en) | A kind of x-ray chest radiograph lung segmentation method and apparatus | |
| Anders et al. | Performance of an atlas-based autosegmentation software for delineation of target volumes for radiotherapy of breast and anorectal cancer | |
| CN111028914B (en) | Artificial intelligence guided dose prediction method and system | |
| US7724936B2 (en) | Image generation apparatus and image generation method for detecting abnormalities | |
| WO2020244172A1 (en) | Plan implementation method and device based on predicted dose guidance and gaussian process optimization | |
| CN103914823B (en) | The method of the quick exact non-linear registration solid medical image based on rarefaction representation | |
| CN104036109A (en) | Image based system and method for case retrieving, sketching and treatment planning | |
| CN108771794B (en) | System and method for generating dose calculations in a radiation treatment plan | |
| CN108721792B (en) | Method, program memory and system for radiation therapy treatment planning | |
| Sjöberg et al. | Multi-atlas based segmentation using probabilistic label fusion with adaptive weighting of image similarity measures | |
| US10512790B2 (en) | Systems and methods for generating radiation treatment plans | |
| CN112601582A (en) | System and method for accelerated on-line adaptive radiotherapy | |
| CN111261296A (en) | Tumor clinical target area automatic delineation method and system based on conditional random vector field | |
| Miura et al. | Impact of deformable image registration accuracy on thoracic images with different regularization weight parameter settings | |
| CN108597589B (en) | Model generation method, target detection method and medical imaging system | |
| CN111888665B (en) | Construction method of three-dimensional dose distribution prediction model based on adaptive countermeasure network | |
| Wu et al. | Development of an accelerated GVF semi-automatic contouring algorithm for radiotherapy treatment planning | |
| CN109564685B (en) | Robust lobe segmentation | |
| CN117992802B (en) | Radiotherapy similarity planning method, system and storage medium based on radiotherapy database |
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 |