WO2012008500A1 - 信号処理装置、信号処理プログラム及び信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体 - Google Patents

信号処理装置、信号処理プログラム及び信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体 Download PDF

Info

Publication number
WO2012008500A1
WO2012008500A1 PCT/JP2011/066007 JP2011066007W WO2012008500A1 WO 2012008500 A1 WO2012008500 A1 WO 2012008500A1 JP 2011066007 W JP2011066007 W JP 2011066007W WO 2012008500 A1 WO2012008500 A1 WO 2012008500A1
Authority
WO
WIPO (PCT)
Prior art keywords
variable vector
target
vector
fluoroscopic image
background
Prior art date
Application number
PCT/JP2011/066007
Other languages
English (en)
French (fr)
Inventor
経康 本間
良尋 高井
春奈 遠藤
慶 市地
正夫 酒井
吉澤 誠
Original Assignee
国立大学法人東北大学
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 国立大学法人東北大学 filed Critical 国立大学法人東北大学
Priority to JP2012524581A priority Critical patent/JP5797197B2/ja
Priority to EP11806830.3A priority patent/EP2604187A4/en
Publication of WO2012008500A1 publication Critical patent/WO2012008500A1/ja
Priority to US13/740,859 priority patent/US8837863B2/en

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/60Editing figures and text; Combining figures or text
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/48Diagnostic techniques
    • A61B6/486Diagnostic techniques involving generating temporal series of image data
    • A61B6/487Diagnostic techniques involving generating temporal series of image data involving fluoroscopy
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5211Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
    • A61B6/5217Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data extracting a diagnostic or physiological parameter from medical diagnostic data
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/54Control of apparatus or devices for radiation diagnosis
    • A61B6/542Control of apparatus or devices for radiation diagnosis involving control of exposure
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/20Analysis of motion
    • G06T7/277Analysis of motion involving stochastic approaches, e.g. using Kalman filters
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5258Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5258Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
    • A61B6/5264Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise due to motion
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61NELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00Radiation therapy
    • A61N5/10X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A61N5/1048Monitoring, verifying, controlling systems and methods
    • A61N5/1049Monitoring, verifying, controlling systems and methods for verifying the position of the patient with respect to the radiation beam
    • A61N2005/1061Monitoring, verifying, controlling systems and methods for verifying the position of the patient with respect to the radiation beam using an x-ray imaging system having a separate imaging source
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10016Video; Image sequence
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10116X-ray image
    • G06T2207/10121Fluoroscopy
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30096Tumor; Lesion

Definitions

  • the present invention relates to a signal processing device, a signal processing program, and a computer-readable recording medium on which the signal processing program is recorded.
  • radiotherapy apparatus that treats an affected area such as a tumor by irradiating radiation such as proton beam, heavy particle beam, or X-ray.
  • radiation such as proton beam, heavy particle beam, or X-ray.
  • the position and shape of the affected area may vary depending on the patient's breathing and the like, and therefore radiation is performed in consideration of the variation.
  • radiation is irradiated only when the affected tissue has moved to a predetermined radiation irradiation position so as to irradiate the affected tissue with radiation while suppressing the influence on the normal tissue of the patient.
  • a method also called follow-up irradiation
  • the affected part is irradiated with radiation intermittently.
  • a method (hereinafter also referred to as tracking irradiation) has been proposed in which the irradiation position of radiation is made to follow the position of a moving affected part.
  • a method has also been proposed in which time fluctuations in the position of the affected area (for example, position fluctuation due to patient breathing) are predicted, and radiation is irradiated to the predicted position.
  • Patent Document 1 discloses a method of acquiring a fluoroscopic image in which a subtle contrast is ensured by controlling the irradiation timing for fluoroscopic image acquisition and the activation timing of therapeutic radiation, or based on time-series data.
  • a method for accurately and safely irradiating the affected part, and the obtained fluoroscopic image of the irradiation target is collated with a reference image for a specific evaluation factor.
  • a method for performing irradiation control with high reliability on the irradiation target is described.
  • a method for grasping (specifying) a change in the position and shape of the affected area there is a method for monitoring the position and shape of the affected area by performing X-ray imaging of the vicinity of the affected area.
  • the X-ray imaging is performed by, for example, an On Board Imager (OBI) mounted in a direction different from the irradiation direction of therapeutic radiation on the radiotherapy device main body, or kV installed independently of the radiotherapy device main body.
  • OBI On Board Imager
  • kV installed independently of the radiotherapy device main body.
  • an X-ray image (so-called fluoroscopic image) lacks resolution and contrast, and the affected tissue and its surrounding tissues are compared to a tomographic image acquired by computed tomography (CT) or the like.
  • CT computed tomography
  • the boundary is unclear.
  • a metal (for example, gold) marker is inserted in advance in the vicinity of the affected area, and X-ray imaging is performed in the vicinity of the affected area including the marker. .
  • variation of an affected part position is estimated by measuring the position of the metallic marker which appears on the fluoroscopic image acquired by X-ray imaging. This is because the metallic marker has sufficient contrast even on a fluoroscopic image.
  • a metal marker when a metal marker is inserted in the vicinity of an affected area, there may be a health risk for the patient.
  • the affected part is a lung tumor or the like and a metal marker is inserted in the vicinity of the lung tumor, emphysema or pneumothorax may occur in about 30% of the entire patient.
  • emphysema or pneumothorax may occur in about 30% of the entire patient.
  • insertion of a metal marker in the vicinity of an affected area is not permitted in the United States and the like. Therefore, it is difficult to measure the position variation of the affected part. As a result, the irradiation accuracy with respect to the affected area is lowered, and a sufficient therapeutic effect cannot be expected, or a side effect may occur due to erroneous irradiation of normal tissue.
  • the present invention has been made in view of the above problems, and an object of the present invention is to measure (estimate) the position, shape, size, and the like of a target portion such as an affected part with high accuracy. Another object is to measure (estimate) the position, shape, size, and the like of a target part such as an affected part with high accuracy without inserting a metallic marker and with a low burden. Furthermore, it is also an object to enable accurate radiation irradiation to the affected area in radiotherapy.
  • the present invention is not limited to the above-described object, and can be positioned as one of other objects that is an effect obtained by each configuration shown in the embodiments to be described later and that cannot be obtained by conventional techniques. .
  • a signal processing device or an image processing device that processes a fluoroscopic image expressed by transparent overlay of a target portion and a background portion different from the target portion, and a plurality of times
  • the fluoroscopic image acquisition unit for acquiring the fluoroscopic image including the target portion, and the fluoroscopic component resulting from the target portion and the background portion resulting from the target portion at a time t of the plurality of times Calculating a component of an image in association with estimated values of distribution information of the target portion and the background portion at the plurality of times, and evaluating the consistency with the fluoroscopic image at the plurality of times;
  • An update unit that updates an estimated value of the distribution information of the target part or the background part at the time t based on the evaluation result of the part. Apparatus can be used.
  • the spatial dimension of the fluoroscopic image is not particularly limited, and may be one-dimensional, two-dimensional, or three-dimensional.
  • the target portion is a detection target in the fluoroscopic image, or a portion corresponding to the target portion in a real space corresponding to the fluoroscopic image.
  • the background portion refers to a portion that is not included in the target portion.
  • One or two or more target portions can be selected. When two or more target portions are selected, the target portions may overlap each other in the fluoroscopic image or may exist apart from each other. Note that there is no essential difference between the target portion and the background portion, and it is possible to regard one or a plurality of target portions as background portions (that is, it is considered that there are a plurality of background portions).
  • the transparent superimposition of the target part and the background part means that the target part and a part or all of the background part are superposed with some weight on each other.
  • the target portion and the background portion cannot be separated only from the fluoroscopic image at a certain time t. Therefore, in the calculation unit, at a certain time t, one or both of the component due to the target portion and the component due to the background portion are distributed to the target portion and the background portion at a plurality of times (which may include time t).
  • the consistency with a fluoroscopic image at a plurality of times (which may be the same as or different from the associated times) is evaluated.
  • the association can be formulated and calculated based on transparent overlay characteristics (for example, recognized as physical characteristics or statistical characteristics).
  • a plurality of time ranges associated with the state at time t can be arbitrarily set. However, if the movement of the target portion and the background portion can be considered to be continuous in time, the state at the time t can be determined if the target portion and the background portion are to be separated with a certain degree of accuracy (depending on the purpose of use). It is preferable to set the plurality of times in a time range in which the change from can be regarded as relatively small (this is referred to as the vicinity of time t).
  • Consistency can be evaluated by, for example, calculating a suitably defined error magnitude quantitatively and comparing it with an appropriate value.
  • the distribution information is information representing the distribution of the target portion or the background portion in the fluoroscopic image.
  • the distribution information can be represented by, for example, spatial position information or luminance information at each position.
  • the position information can be, for example, presence information in each region (which may be composed of one pixel or plural pixels) of interest as a processing target, and luminance information. Can be information on luminance values representing each region.
  • the position information is handled as linear, it can be expressed as a coefficient matrix having coefficients for pixel positions in each region, for example.
  • the target portion and the background portion may be separated by performing a statistical analysis.
  • the validity of separation is calculated by calculating the presence / absence / degree of statistical independence between the target part and the background part for at least one of the position, shape, and brightness of the target part and the background part. Can be evaluated.
  • An example of such an analysis is an independent component analysis.
  • a processing apparatus which is a fluoroscopic image acquisition unit that acquires a fluoroscopic image for the target portion at a certain time t, and M ⁇ N (M and N are natural numbers) of the fluoroscopic images acquired by the fluoroscopic image acquisition unit.
  • K fluoroscopic images acquired at the time t and K times
  • K is a natural number
  • t 1 , t 2 ,..., T K before the time t.
  • Observation that generates an observation vector b (t, t 1 , t 2 ,..., T K ) for (K + 1) time which is composed of luminance values of (K + 1) ⁇ M ⁇ N regions of interest.
  • a vector generator and a variable vector composed of luminance values attributed to each of the L target portions in the region of interest.
  • a variable vector update unit that updates the target variable vector I a (t p ) and the background variable vector I b (t p ), respectively, the observation vector b, and the target variable vector updated by the variable vector update unit
  • a variable vector I composed of I a and a background variable vector I b , positions of the L target portions and background portions at the time t estimated using the variable vector I, and the K past times t 1 , t 2 ,..., t K were similarly estimated From a coefficient matrix A of (K + 1) MN ⁇ (L + 1) MN relating to the positions of the L target parts and the background part at each time, f is a deterministic or stochastic function. ,formula
  • variable vector calculation unit that calculates the variable vector I so as to enhance the evaluation of the evaluation value PI defined by the acquisition of the fluoroscopic image in the fluoroscopic image acquisition unit, and the observation vector generation unit Signal processing configured to generate the observation vector b, update the variable vector I in the variable vector update unit, and calculate the variable vector I in the variable vector calculation unit while shifting the time
  • An apparatus or an image processing apparatus can be used.
  • a fluoroscopic image acquisition function for acquiring the fluoroscopic image including: a component of the fluoroscopic image due to the target portion and a component of the transparent image due to the background portion at a certain time t among the plurality of times; Calculated in association with estimated values of the distribution information of the target portion and the background portion at the plurality of times, and a calculation function for evaluating the consistency with the fluoroscopic images at the plurality of times, and an evaluation result by the calculation function
  • a computer-readable recording medium in which the program is recorded can be used.
  • a signal processing method having steps for performing each of the above processes can be used.
  • the present invention it is possible to measure (estimate) the position, shape, size, and the like of a target portion such as an affected part with high accuracy.
  • a target portion such as an affected part
  • radiation therapy it is possible to accurately irradiate the affected area.
  • FIG. It is a figure which shows an example of a radiotherapy apparatus. It is a functional block diagram of the control part shown in FIG. It is a figure which shows an example of an optical flow. It is a figure which shows an example of a block matching (BM) method. It is a figure which shows the difference
  • FIG. It is a figure which shows an example of the signal processing method of this example. It is a figure which shows an example of the signal processing method of this example.
  • (A) is a figure which shows an example of the correct image of object (tumor) part
  • (B) is a figure which shows an example of the correct image of a background part
  • (C) shows an example of the fluoroscopic image of phantom data.
  • FIG. It is a figure which shows an example of the fluoroscopic image of phantom data. It is a figure which shows the object part isolate
  • (A) to (C) are diagrams comparing the target image separated from the fluoroscopic image and the background image. It is a figure which shows an example of the fluoroscopic image of case data. It is a figure which shows the object part isolate
  • FIG. 1 illustrates a radiotherapy device. 1 includes, for example, a control unit 1, a radiation irradiation device 90, a measurement radiation irradiation unit 305, and a sensor unit 306.
  • a control unit 1 controls the radiation therapy device
  • the affected part 260 can be irradiated with therapeutic radiation.
  • the position, shape, size, and the like of the affected part 260 may fluctuate with time due to the patient's respiratory motion and the like. Therefore, in this example, the measurement radiation irradiation unit 305 irradiates the periphery of the affected part 260 with radiation for photographing a fluoroscopic image (for example, X-rays), and the sensor part 306 transmits the radiation for photographing the fluoroscopic image that has passed through the affected part 260. By detecting, a radioscopic image around the affected part 260 (hereinafter, also simply referred to as a fluoroscopic image) is taken and input to the control unit 1.
  • a fluoroscopic image for example, X-rays
  • the control unit 1 measures (estimates) temporal variations such as the position, shape, and size of the target portion such as the affected part 260 in the fluoroscopic image. Moreover, the control part 1 may control the irradiation position and range of the therapeutic radiation irradiated from the radiation irradiation apparatus 90 based on the said measurement result. At this time, the therapeutic radiation irradiated from the radiation irradiation device 90 is used as a substitute for measurement radiation, and a fluoroscopic image of the affected area 260 is taken and measured using an electronic portal imaging device (EPID) or the like using the radiation that has passed through the affected area 260 ( Estimation). Alternatively, both may be photographed and measured (estimated). Further, the control unit 1 may control the bed 230 to match the affected part 260 with the radiation irradiation position.
  • EPID electronic portal imaging device
  • the plurality of fluoroscopic images may be taken by the same kind of fluoroscopic image device or may be taken by different types of devices.
  • the 3D measurement with movement may be performed by fusing the 4D-CT data obtained by photographing the movement of the affected part 260 in advance and the fluoroscopic image data.
  • the position, shape, size, and the like of the affected part 260 were basically measured using a fluoroscopic image acquired by an X-ray simulator.
  • the X-ray simulator can take a fluoroscopic image of 2 frames per second.
  • the On Board Imager can take a fluoroscopic image of 15 frames per second, enabling more accurate measurement. Can be done.
  • control unit 1 includes, for example, a signal processing unit 2 and a radiation irradiation control unit 3 as shown in FIG.
  • the signal processing unit (signal processing device) 2 measures (estimates) the position of the affected part 260 that varies with time based on the fluoroscopic images taken by the measurement radiation irradiation unit 305 and the sensor unit 306.
  • the signal processing unit 2 includes, for example, a fluoroscopic image acquisition unit 4, an observation vector generation unit 5, a variable vector update unit 6, a variable vector calculation unit 7, and an estimation unit 8.
  • the operations of the fluoroscopic image acquisition unit 4, the observation vector generation unit 5, the variable vector update unit 6, the variable vector calculation unit 7 and the estimation unit 8 will be described in detail later in (1.6). An outline of the operation of each unit will be described.
  • the fluoroscopic image acquisition unit 4 acquires a fluoroscopic image photographed by the measurement radiation irradiation unit 305 and the sensor unit 306.
  • the fluoroscopic images are acquired at a plurality of times including a certain time t, and the fluoroscopic images generally include L (L is a natural number) target portions that are measurement (estimation) targets such as position, shape, and size.
  • the fluoroscopic image can be considered as (L + 1) objects because there are L target parts and a background part other than the target part (however, there is no essential difference between the target and the background).
  • it can be said that it is expressed by transparent superposition with a plurality of backgrounds in which some objects are considered as backgrounds.
  • the observation vector generation unit 5 has M ⁇ N (M and N are natural numbers) pixels (pixels digitized in an analog image) or a plurality of pixels of the fluoroscopic image at the time t acquired by the fluoroscopic image acquisition unit 4.
  • a fluoroscopic image is stored by paying attention to a certain area consisting of some representative value, and the brightness value (hereinafter also referred to as density) of the MN target areas and the past K from the time t (K is a natural number) )
  • Pieces of K perspective images acquired and stored in the same manner as described above at times t 1 , t 2 ,..., T K , together with an observation vector b (t consisting of (K + 1) MN regions of interest. , t 1 , t 2 ,..., t K ).
  • the fluoroscopic image may be stored in the fluoroscopic image acquisition unit 4.
  • the variable vector updating unit (updating unit) 6 includes variable vectors I a1 , I a2 ,..., I aL composed of luminance values caused by L target portions in the region generated by the observation vector generating unit 5.
  • L number of variables and variable vector I b where a luminance value having the background portion in the generated area, related to the interest region of the observation vector b and the time t fluoroscopic image even at a certain time t p of the past from Update using the target variable vector I a (t p ) and the background variable vector I b (t p ) consisting of the vectors I a1 (t p ), I a2 (t p ), ..., I aL (t p ), respectively To do.
  • the variable vector calculation unit (calculation unit) 7 includes the observation vector b generated by the observation vector generation unit 5, the target variable vector I a and the background variable vector I updated by the variable vector update unit 6. b, the position of the L target parts and the background part at the time t estimated using the variable vector I, and the K past times t 1 , t 2 ,. from the L target portion and a position of the background portion (K + 1) MN ⁇ ( L + 1) MN of the coefficient matrix and (function) a at each time estimated similarly in K, determine f As a theoretical or probabilistic function
  • the variable vector I is calculated so that the evaluation value PI defined by (1) is an appropriate value such as an optimum value.
  • the estimation unit 8 calculates the time of the target portion based on the target variable vector I a set with an initial value by the variable vector update unit 6 and the target variable vector I a updated by the variable vector update unit 6.
  • the target variable vector I a before being updated by the variable vector update unit 6 and the target variable vector I a after being updated by the variable vector update unit 6 may be estimated. Based on the above, the temporal position variation of the target portion may be estimated. Further, based on the updated target variable vector I a, may be estimated on the shape and size of the subject.
  • the signal processing unit 2 acquires the fluoroscopic image in the fluoroscopic image acquisition unit 4, generates the observation vector b in the observation vector generation unit 5, and the variable vector I in the variable vector update unit 6. And the calculation of the variable vector I by the variable vector calculation unit 7 are performed while shifting the time. Information such as the position, shape and size of the target part such as the affected part 260 measured by the signal processing unit 2 is output to the radiation irradiation control unit 3.
  • the radiation irradiation control unit 3 controls the irradiation position and irradiation range of the therapeutic radiation irradiated from the radiation irradiation apparatus 90 or the position of the bed 230 based on the information from the signal processing unit 2.
  • a general position estimation method was tried for the above-described fluoroscopic image. This is to extract a movement amount vector (u, v) for each pixel of a target portion (for example, an affected part 260) T between each frame (also referred to as a moving image frame) of a fluoroscopic image taken at different times. is there.
  • FIG. 3 is an example in which the movement amount vector (u, v) is obtained for the target portion T at time t (> 0) and the target portion T at a discrete next time (t + 1).
  • the next time (t + 1) may be an arbitrary time in the future than t. However, for the sake of simplicity, the next time (t + 1) is simply (t + 1).
  • BM general block matching
  • a region of a certain size including a target point of a target portion is used as a template block, and the next frame (at (t + 1)) captured at time (t + 1).
  • the evaluation function such as the difference with the template block of the previous frame is optimized. This is a method of extracting the movement displacement of each point constituting the target portion between frames by using the points having the same relative position in the block as corresponding points.
  • a movement amount vector (u, v) is calculated.
  • the luminance value of a certain pixel at the position (m, n) of the target portion in the frame photographed at time t is I (m, n, t), and the position of the extraction target in the frame photographed at time (t + 1)
  • the luminance value of the pixel at (m, n) is I (m, n, t + 1)
  • the luminance value I (m, n, t) of a certain pixel in the previous frame and the position of the pixel in the next frame The least square error Mean Square Error (MSE) with the luminance value I (m + u, n + v, t + 1) of the pixel moved by the movement amount vector (u, v) is defined by the following equation (a).
  • ⁇ m, n ⁇ is a set of all pixel positions existing in the M ⁇ N block.
  • a flow (U, V) that minimizes this MSE is calculated for each block by the following equation (b).
  • the difference (estimated error) between the actual position of the affected area 260 and the estimated position was about 1.2 [ mm] to about 2.4 [mm].
  • a contour having high contrast by visual observation among the tumor part of the target affected part 260 The part was set manually as a point of interest.
  • the block size was square N ⁇ N. In FIG. 5, the horizontal axis represents the block size (N).
  • the BM method when used for estimating the position of the affected part 260, the accuracy required clinically (for example, the estimation error is within 1 [mm]) in spite of the improvement in accuracy as described above. Cannot be achieved. (1.4) Problems of general position estimation method As described above, the BM method cannot achieve a desired estimation accuracy (error is within 1 [mm]). Two are possible.
  • the fluoroscopic images are very unclear, and the luminance value may be different between frames due to noise or the like even for the same target.
  • the luminance value of the affected part 260 is represented by an integrated value of the luminance value of the target part itself such as the affected part 260 and the luminance value of the background part. That is, the luminance value of the target portion does not change even when the target portion moves. Since this condition is also premised on other general methods other than optical flow extraction such as phase-only correlation method, the failure to satisfy the condition is a more essential cause of a decrease in accuracy.
  • Template matching In the BM method, as illustrated in FIG. 6, at least one of the previous frame and the next frame is used to calculate an evaluation value such as a difference in the target portion between the previous frame and the next frame.
  • an evaluation value such as a difference in the target portion between the previous frame and the next frame.
  • a template matching (TM) method can be used instead of the BM method.
  • the TM method is a method in which a template of a target portion is acquired in advance and a range having the highest correlation with the template is searched for in each frame at each time.
  • a template used in the TM method can be acquired by, for example, defining a range having a luminance value equal to or higher than a predetermined value in a frame at a certain time as a template.
  • the template may be determined based on the position, shape, size, and the like of the affected part 260 measured by another measuring device (CT, MRI (magnetic resonance imaging, etc.)).
  • the least square error (MSE) between the luminance value I temp (m, n) of the template and the luminance value I (m + u, n + v, t) of the target portion in the frame at time t is expressed by the following equation (c ).
  • the position of the target portion is measured (estimated) based on the evaluation value of the target portion template that is not affected by noise and the target portion in each frame.
  • the fluoroscopic image is less susceptible to ambiguity.
  • the density distribution over the details of the target portion affects the evaluation value.
  • the TM method for example, by using a template having a constant density distribution, the difference in the details of the target portion is not significantly affected. Can be expected.
  • it is also effective to add a noise term to equation (c) and perform statistical processing.
  • the apparent luminance value of the target part such as the affected part 260 is the original luminance value of the target part and the affected part 260. It is determined by the integrated value with the luminance value of the background part such as the peripheral organs of Therefore, even if attention is paid to a part of the target part, if the luminance value of the background part overlapping the part is temporally different, the apparent luminance value of the part is not constant over time. Therefore, generally, the preconditions of the conventional optical flow method described above are not satisfied in the fluoroscopic image.
  • the luminance value of the target portion does not change according to the luminance value of the background portion. That is, when attention is paid to a part of the target part in the reflected image, the apparent luminance value of the part is constant even if the luminance value of the background part overlapping the part is temporally different. Therefore, generally, in the reflection image, the preconditions of the above-described conventional optical flow method are satisfied.
  • the apparent (observed) luminance value I (m, n, t) in the fluoroscopic image is converted into the original luminance value I a of the target portion that moves in time. (M, n, t) and a luminance value I b (m, n, t) of a background part that moves differently from the target.
  • the background displacement is set to 0 (that is, still)
  • the luminance value I observed in the fluoroscopic image is By estimating the position of the original luminance value I a (m, n, t) that is not temporally moving (m, n, t) but using the TM method or the like, the target portion is similar to the reflected image. Can be measured (estimated) with high accuracy.
  • the separation accuracy of Ia using the estimated position information not only the position of the target part but also the shape and size are measured (estimated) by repeating each process. ) Is possible.
  • a luminance value vector (hereinafter also referred to as an observation vector) I (m, n) of a certain region (also referred to as a region of interest (ROI)) of a fluoroscopic image at time t. , t) is the luminance value vector I a (m of the target portion, n, t) and the luminance value vector I b (m background section, n) and by, can be expressed by the following equation (e) .
  • the equation represented by the equation (e) is indefinite in which the number of equations independent of the number of unknown variables is insufficient, and the variable vector I a (I a22 to I a33 ) of the target portion and The distribution (I b22 to I b33 ) of the variable vector I b in the background portion is not uniquely determined. This means that an element part constituting the density cannot generally be separated from one fluoroscopic image.
  • the observation frames acquired at time t and the observation frames acquired at another time are temporally integrated, and each ROI the observation vector I, the variable vector I a of the target portion by the variable vector I b of the background portion, expressed by the following equation (f).
  • the matrix on the left side of the left side is a coefficient matrix (function) A, for example, a matrix determined by the position of the target portion at different times t and (t + 1).
  • the matrix on the right side of the left side is a variable vector I indicating the density of the ROI target part and the background part at time t.
  • the matrix on the right side is an observation vector b indicating the luminance value observed in the ROI of the fluoroscopic image acquired at different times t and (t + 1).
  • the fluoroscopic image is composed of 4 ⁇ 4 16 pixels or 16 representative values of a plurality of pixels
  • the variable vector I is merely an 8 ⁇ 1 matrix, and the number of pixels (or representative values) of the fluoroscopic image, the number of pixels (or representative values) of the ROI, the number of target portions, the past The size of each matrix changes according to the number of times.
  • the coefficient matrix A is full rank, the variable vector I associated with the coefficient matrix A is uniquely obtained.
  • the coefficient matrix A in particular, the coefficient matrix A 5, 6, 7, 8.
  • the coefficient matrix A is estimated.
  • each element of the variable vector I associated with the coefficient matrix A can be set, for example, as the initial value of the luminance value of the template created by the template creation method described above, but is only an assumed value. Generally different from the actual value. Therefore, in this example, the above-described coefficient matrix A and variable vector I are started from a certain initial value, and updating using mutual estimated values is repeated to improve the estimation accuracy.
  • FIG. 12 is an example of a flowchart showing the operation of this example.
  • the signal processing unit 2 assigns “0” to the time parameter t as an initial setting, and sets a region of interest (ROI) in the fluoroscopic image.
  • the variable vector update unit 6 sets, for example, the luminance value vector I a template of the target portion template as the initial value of the variable vector I a (step A1).
  • the ROI is set so as to include a target part such as the affected part 260.
  • a target part such as the affected part 260.
  • an area that always includes the target portion is acquired automatically, manually by visual observation, or in advance.
  • the variable vector updating unit 6 roughly outlines the contour of the target portion, for example, automatically based on a CT or MRI image obtained by separately capturing the target portion, or manually by visual inspection. Assuming that a portion (range) having a luminance value greater than or equal to a predetermined value is a target portion in a fluoroscopic image extracted or acquired at a certain time, the portion can be set as a template for the target portion. . Furthermore, I a template can be calculated by setting the luminance value of the extracted portion to an appropriate constant value, or based on information about the brightness of the template .
  • the fluoroscopic image acquisition unit 4 acquires a fluoroscopic image for the target portion at time t. That is, the fluoroscopy image acquisition unit 4 and the observation vector generation unit 5 acquire the brightness value vector I (t) in the ROI set in step A1 for the fluoroscopic image acquired at time t (step A2).
  • the variable vector update unit 6 is an estimated value of the variable vector I that is used to estimate the movement amount vector (u, v) of the target portion based on the following steps A3 to A8.
  • the variable vector update unit 6 initializes the variable vector I b ′ (Step A4).
  • the observation vector generation unit 5 may measure the brightness per unit area in each block in the ROI of the fluoroscopic image, and calculate the luminance value based on the measurement result. Then, the variable vector update unit 6 updates I b ′ (step A7). I b at step A7 'update, for example, before time variable vector calculated by I b by the variable vector calculation unit 7' that the value of is set to variable vector I b 'at the next time Done.
  • variable vector update unit 6 updates the variable vector I a ′ (step A8).
  • the update of I a ′ in step A8 is performed by, for example, the variable vector I b ′ initialized or updated in step A4 or A7 from the observation vector I (t) generated by the observation vector generation unit 5 at time t. Is performed, and the subtraction result is set to I a ′.
  • the estimation unit 8 applies, for example, the TM method or the gradient method between I a template set in step A1 and the variable vector I a ′ updated in step A8, thereby moving the movement amount of the target portion.
  • Vector estimated values (u ′ (t), v ′ (t)) are calculated (estimated) (step A9).
  • the target part is a rigid body in which all the points constituting the target part have the same movement amount, or all or some of the points are other points. May be strictly estimated as a non-rigid body having a different movement amount. For the sake of simplicity, the following description will be made assuming that the object is a rigid body.
  • the estimation unit 8 for example, by applying the BM method or gradient method between the previous variable vector I a and 'a variable vector I a that has been updated in step A8' to be updated in step A8
  • the estimated value (u ′ (t), v ′ (t)) of the movement vector of the target portion may be calculated (estimated), or the variable vector I a after being updated in I a template or step A8. It may be calculated (estimated) by a phase-only correlation method with 'or a method using a particle filter.
  • variable vector calculation unit 7 estimates and updates the coefficient matrix A based on the estimated value (u ′ (t), v ′ (t)) of the movement amount vector estimated by the estimation unit 8. That is, the variable vector calculation section 7 (function updating unit), based on the variable initial values are updated with the set variable vector I a vector I a, a function of estimating and updating the coefficient matrix (Function) A Have. Alternatively, the variable vector calculation unit 7, based on the variable vector I a after being the update to the previous variable vector I a to be updated has a function of estimating and updating the coefficient matrix A.
  • variable vector calculation unit 7 calculates I (t) by calculation so as to minimize the square of each element of the error e defined by the following equation (g), for example (step A10).
  • the separation concentration is estimated using an optimization method or the like instead of strictly obtaining a solution.
  • the method of minimizing the square of the error e includes so-called computational intelligence such as an optimization method such as steepest descent method, conjugate gradient method, Newton method, artificial neural network, fuzzy theory, evolutionary calculation theory, etc. Techniques can also be used. It is also possible to use a stochastic optimization method by adding a statistical noise term to the error function. This is an effective method for an unclear image having a lot of noise. Furthermore, since there are many images resulting from the tumor part that is the target in this example and only the target exists and the brightness other than the target part is 0, the brightness value is greater than a certain value for the estimated I ′ a Alternatively, noise may be removed from a small pixel value by threshold processing or the like. This can be expected to reduce the accumulation of estimation errors especially during the optimization process.
  • computational intelligence such as an optimization method such as steepest descent method, conjugate gradient method, Newton method, artificial neural network, fuzzy theory, evolutionary calculation theory, etc. Techniques can also be used. It is also possible to use a stochastic optimization method by adding a statistical noise term
  • variable vector calculation unit 7 calculates the variable vector I so as to make any evaluation index such as the error e appropriate by using an optimization method based on the estimated coefficient matrix A.
  • the signal processing unit 2 increments the value of time t (step A11), and repeatedly executes the processing of step A2 to step A11. That is, the signal processing unit 2 in this example acquires the fluoroscopic image in the fluoroscopic image acquisition unit 4, generates the observation vector b in the observation vector generation unit 5, and the variable vector I a in the variable vector update unit 6. and updating of I b, and the calculation of the variable vector I a and I b in variable vector calculation unit 7 configured to perform while shifting the time.
  • variable vector update unit 6 updates the variable vectors I a and I b using the variable vectors I a and I b calculated by the variable vector calculation unit 7, respectively.
  • the target portion and the background portion can be separated from the fluoroscopic image, and the position of the target portion is increased by estimating the position of the separated target portion. It becomes possible to measure (estimate) accurately.
  • the shape and size of the target portion can be measured (estimated) with high accuracy by separation.
  • the phantom data is a transparent image as shown in FIG. 13B with no target tumor having a ROI size of 100 ⁇ 100, and an artificially moving tumor image as shown in FIG.
  • FIG. 14 shows a change example of the phantom data frame with time development.
  • 12 frames of fluoroscopic images appropriately selected from a total of 110 frames are arranged from upper left to lower right in chronological order.
  • the tumor moves mainly in the vertical direction around the center of the ROI, simulating the movement of an actual lung tumor.
  • the signal processing method of this example is applied to each frame shown in FIG. 14, even if an initial value considerably different from the correct answer is manually set appropriately as shown in FIGS. In both cases, the target part and the background part could be separated accurately.
  • FIG. 15 it can be seen that the tumor image shown in FIG. 13 (A) can be separated and extracted from the initial value appropriately manually set by visual inspection with time development.
  • the initial value error was 300 pixels or more out of 10,000 ROI sizes and the density difference was 20 per pixel, whereas in the final frame, the error was 0 pixel and the density difference per pixel was about 0.1 on average. It was successfully reduced to the following, indicating that separation was possible with high accuracy.
  • FIG. 17 shows an example in which a separation result obtained by repeating initial value setting and updating is compared with a correct image.
  • the initial image (B) the shape and size of the tumor part is significantly different from the correct image (A), and the effect is clearly seen in the background image, but in the separation result image (C), the shape and size of the tumor are also accurate. Can be separated. Although some errors remain in the density distribution, further accurate estimation can be performed by repeatedly applying the signal processing method of this example. It is difficult to see the effect.
  • the estimation error is almost zero, and an accurate position estimation result is obtained. It was.
  • the conventional BM method was applied to the same phantom data, the average error was about 3 [mm] due to the high contrast rib contour of the background image. Can be confirmed.
  • FIG. 18 shows an example of a fluoroscopic image before separating the target portion and the background portion of ROI size 60 ⁇ 60.
  • continuous perspective images from the first frame to the twelfth frame are arranged from the upper left to the lower right in time series order.
  • the signal processing method of this example was applied to each frame shown in FIG. 18, the target portion and the background portion could be separated as shown in FIGS.
  • the estimation error is 0.61 ⁇ 0.18 [mm]. It became. Moreover, the shape and size of the target part could be measured with high accuracy. Although the phantom data was updated for 100 frames, the case data used was limited to use only 12 frames, so the accuracy was somewhat lower than that of the phantom data, but clinically sufficient high-accuracy measurement was possible. It has become. In actual clinical practice, it is possible to acquire data of about 100 frames by performing observations for about 1 or 2 minutes before treatment irradiation, so even more accurate separation and position measurement is possible for case data. Is possible.
  • the estimation error was It was 0.81 ⁇ 0.04 [mm].
  • the accuracy is improved as compared with the BM method, and the effect of the TM method can be confirmed.
  • the TM method in this case can be expected to have high accuracy if the contrast is high, such as the outline of the target portion being clear as in the BM method, it is a trial-and-error approach with the contour portion having a high contrast visually as a point of interest.
  • the ROI size with the smallest estimation error was adopted. Since manual setting by visual inspection is a cumbersome task requiring great care, performing it in a busy clinical setting places a heavy burden on the setter. Therefore, it is difficult to perform highly accurate measurement without separating the target and the background.
  • the target portion and the background portion are separated from the fluoroscopic image, and the position of the separated target portion is estimated, so that the target portion is obtained from the fluoroscopic images acquired at different times. It is possible to measure (estimate) the position, shape, size, and the like of the object with high accuracy.
  • the position, shape, size, etc. of the target portion can be measured (estimated) with high accuracy and safely without inserting a metallic marker into the observation target. Can do.
  • I ′ can be calculated by defining the error e by the following equation (i) and applying an optimization method such as the steepest descent method to the error e (see the circled number 3 in FIG. 21). .
  • a radiation generator pulse modulator 250, klystron 100, cooling device 110, waveguide 120, vacuum pump 130, etc.
  • a radiation generator pulse modulator 250, klystron 100, cooling device 110, waveguide 120, vacuum pump 130, etc.
  • the signal processing device signal processing unit 2
  • the irradiation position and irradiation range of the radiation are calculated, and the calculation result is obtained.
  • a radiotherapy apparatus comprising a drive control unit (control unit 1) for driving and controlling the collimator unit based on Erareru.
  • the function as the signal processing unit 2 described above may be realized by a computer (including a CPU, an information processing device, and various terminals) executing a predetermined application program (signal processing program). That is, in the signal processing apparatus for estimating the temporal position variation of the target part in the fluoroscopic image expressed by transparent overlay of the target part and a background part different from the target part, A signal processing program for causing a computer to realize an estimation function, wherein the signal processing program is acquired by a fluoroscopic image acquisition function for acquiring a fluoroscopic image for the target portion at a certain time t and the fluoroscopic image acquisition function.
  • variable vector I consisting of and before estimated using the variable vector I Past time t 1 position and of said K of the L pieces of object portion and background portion at time t, t 2, ..., the at respective times estimated similarly in t K L pieces of object portion and background From the (K + 1) MN ⁇ (L + 1) MN coefficient matrix A with respect to the position of the part, let f be a deterministic or stochastic function,
  • variable vector calculation function for calculating the variable vector I so as to make the evaluation value PI defined by, for example, an appropriate value, for example, obtaining the fluoroscopic image by the fluoroscopic image acquisition function,
  • the generation of the observation vector b by the observation vector generation function, the update of the variable vector I by the variable vector update function, and the calculation of the variable vector I by the variable vector calculation function are performed while shifting the time. It is an example of the signal processing program which operates a computer.
  • the above program is, for example, a flexible disk, CD (CD-ROM, CD-R, CD-RW, etc.), DVD (DVD-ROM, DVD-RAM, DVD-R, DVD-RW, DVD + R, DVD + RW, etc.). Or the like recorded in a non-transitory computer-readable recording medium.
  • the computer can read the signal processing program from the recording medium, transfer it to the internal storage device or the external storage device, store it, and use it.
  • the program may be recorded in a storage device (recording medium) such as a magnetic disk, an optical disk, or a magneto-optical disk and provided from the storage device to a computer via a communication line.
  • the computer is a concept including hardware and an OS (operating system), and means hardware operating under the control of the OS. Further, when the OS is unnecessary and the hardware is operated by the application program alone, the hardware itself corresponds to the computer.
  • the hardware includes at least a microprocessor such as a CPU and means for reading a computer program recorded on a recording medium.
  • the application program as the signal processing program includes program code for causing the computer as described above to realize the function as the signal processing unit 2. Also, some of the functions may be realized by the OS instead of the application program.
  • the recording medium according to this embodiment includes an IC card, ROM cartridge, magnetic tape, punch card, and internal storage device of a computer ( It is also possible to use various computer-readable media such as a memory such as a RAM and a ROM, an external storage device, and a printed matter on which a code such as a barcode is printed.
  • each structure and each process of the radiotherapy apparatus mentioned above and the signal processing part 2 may be selected as needed, and may be combined suitably.
  • the position estimation method when the target portion has moved out of the ROI has been described.
  • the solution cannot be calculated by inverse matrix calculation, and a method using an optimization method while correcting an erroneous constraint condition shown in the modification is used.
  • the signal processing method of the present example is, for example, an estimation of the position, shape, size, etc. of a blood vessel imaged by angiography, and a movement of a part of or all of a plurality of permeable microorganisms.
  • the present invention can be applied to non-destructive inspection.
  • the position of the affected part 260 may be predicted using information on the position, shape, size, etc. of the target part such as the affected part 260 estimated by the signal processing method of this example as an input signal.

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Physics & Mathematics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Optics & Photonics (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Multimedia (AREA)
  • Physiology (AREA)
  • Image Analysis (AREA)
  • Radiation-Therapy Devices (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

 患部等の対象部分と前記対象部分とは異なる背景部分との透過的な重ね合わせにより表現される透視画像を処理する信号処理装置において、複数の時刻について、前記対象部分を含む前記透視画像を取得する透視画像取得部と、前記複数の時刻のうちのある時刻tにおける前記対象部分に起因する成分及び前記背景部分に起因する透過画像の成分を、前記複数の時刻における前記対象部分及び前記背景部分の分布情報の推定値と関連づけて計算し、前記複数の時刻における前記透視画像との整合性を評価する計算部と、前記計算部の評価結果に基づいて、前記時刻tにおける前記対象部分又は前記背景部分の分布情報の推定値を更新する更新部と、を備える。

Description

信号処理装置、信号処理プログラム及び信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体
 本発明は、信号処理装置、信号処理プログラム及び信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体に関する。
 例えば、腫瘍などの患部に、陽子線や重粒子線またはX線等の放射線を照射して治療を行なう放射線治療装置がある。
 上記の放射線治療装置により患部に放射線を照射する場合、患者の呼吸などにより患部の位置や形状が変動することがあるため、前記変動を考慮して放射線照射を行なう。
 このような放射線治療の一例として、例えば、患者の正常組織への影響を抑制しつつ患部組織へ放射線を照射すべく、患部組織が放射線の所定の照射位置に移動した場合にだけ放射線を照射する方法(追撃照射とも称する)が実用化されている。つまり、追撃照射では、患部に対して断続的に放射線が照射される。 
 また、放射線の照射位置を、移動する患部位置に追随させて照射する方法(追尾照射とも称する)が提案されている。さらに、患部位置の時間変動(例えば、患者の呼吸による位置変動)を予測して、予想された位置に放射線を照射する方法も提案されている。
 例えば、下記特許文献1には、透視画像取得用の照射タイミングおよび治療用放射線の起動タイミングを制御することにより、微妙な濃淡差が確保された透視画像を取得する方法や、時系列データに基づいて治療放射線照射時点での最も確からしい患部位置を予測することにより、患部に精度良く安全に放射線を照射する方法、また、取得された照射対象の透視画像を特定の評価因子について基準画像と照合することにより、照射対象に対して信頼性の高い照射制御を行なう方法などが記載されている。
特開2006-51199号公報
 患部の位置や形状の変動を把握(特定)する方法として、患部近傍をX線撮影することにより、患部の位置や形状をモニタリングする方法がある。なお、前記X線撮影は、例えば、放射線治療装置本体に治療用放射線の照射方向とは異なる方向に搭載されるOn Board Imager(OBI)や、放射線治療装置本体とは独立して設置されるkV-X線透視装置、あるいは治療用放射線を併用するMV-X線透視装置などにより行なわれる。
 しかしながら、一般的に、X線撮影された画像(いわゆる透視画像)では、解像度やコントラストが不足するほか、Computed Tomography(CT)などにより取得される断層像に比して、患部組織とその周辺組織との境界が不明瞭である。
 このため、従来の方法では、患部位置を容易に把握すべく、患部近傍に金属製(例えば金など)のマーカを予め刺入しておき、当該マーカを含んだ患部近傍のX線撮影を行なう。そして、X線撮影により取得される透視画像上に表れる金属性マーカの位置を計測することにより、患部位置の変動を推定する。これは、金属性マーカが透視画像上でも十分なコントラストを有するからである。
 しかしながら、患部近傍に金属製マーカを刺入する場合、患者に健康上のリスクが発生することがある。例えば、患部が肺腫瘍などであり、肺腫瘍の近傍に金属製マーカが刺入される場合、患者全体の約30%に肺気腫や気胸などが発生することがある。
 また、上記リスクの観点から、米国などでは患部近傍への金属製マーカの刺入は認められていない。そのため、患部の位置変動を計測することが困難である。結果、患部に対する放射線の照射精度が低下し、十分な治療効果が期待できなかったり、正常組織への誤照射により副作用が発生したりすることがある。
 本発明は、以上のような課題に鑑みたもので、患部などの対象部分の位置,形状及び大きさなどを高精度に計測(推定)することを目的とする。
 また、金属性マーカを刺入することなく安全かつ低負担で、患部などの対象部分の位置,形状及び大きさなどを高精度に計測(推定)することも目的の一つである。
 さらに、放射線治療において、患部への正確な放射線照射を可能とすることも目的の一つである。
 なお、前記目的に限らず、後述する実施形態に示す各構成により導かれる作用効果であって、従来の技術によっては得られない作用効果を奏することも他の目的の一つとして位置付けることができる。
 (1)第1の案として、対象部分と前記対象部分とは異なる背景部分との透過的な重ね合わせにより表現される透視画像を処理する信号処理装置または画像処理装置であって、複数の時刻について、前記対象部分を含む前記透視画像を取得する透視画像取得部と、前記複数の時刻のうちのある時刻tにおける前記対象部分に起因する前記透視画像の成分及び前記背景部分に起因する前記透視画像の成分を、前記複数の時刻における前記対象部分及び前記背景部分の分布情報の推定値と関連づけて計算し、前記複数の時刻における前記透視画像との整合性を評価する計算部と、前記計算部の評価結果に基づいて、前記時刻tにおける前記対象部分又は前記背景部分の分布情報の推定値を更新する更新部とを備える、信号処理装置または画像処理装置を用いることができる。
 ここで、透視画像の空間次元は特に限定されるものではなく、1次元、2次元、3次元のいずれであってもよい。また、対象部分とは、透視画像内において検出対象となるもの、または、透視画像に対応する現実空間等において対象部分に対応するものをいう。さらに、背景部分は、対象部分に含まれない部分をいう。また、対象部分は、一又は二以上を選ぶことが可能であり、二以上を選ぶ場合において、対象部分は、透視画像内において重ね合わさっていてもよいし、離れて存在していてもよい。なお、対象部分と背景部分とには本質的な違いはなく、一または複数の対象部分を背景部分とみなす(つまり、複数の背景部分が存在すると考える)ことも可能である。
 また、対象部分と背景部分との透過的な重ね合わせとは、対象部分と背景部分の一部又は全部とが、互いになんらかの重みをもって重ね合わされることをいう。透過的な重ね合わせが行なわれる場合、ある時刻tの透視画像だけからは、対象部分と背景部分とを分離することはできない。そこで、計算部では、ある時刻tにおいて、対象部分に起因する成分と背景部分に起因する成分の一方または両方とを、複数の時刻(時刻tを含んでもよい)における対象部分及び背景部分の分布情報と関連づけ、複数の時刻(関連づけた複数の時刻と同じであっても、違っていてもよい。)の透視画像との整合性を評価する。ここで、関連づけは、透過的な重ね合わせの特性(例えば、物理的特性や統計的特性などとして認識される)に基づいて定式化し演算することができる。時刻tの状態と関連付けられる複数の時刻の範囲は、任意に設定可能である。しかし、対象部分及び背景部分の動きが、時間的に連続であるとみなせる場合において、対象部分と背景部分とを(利用目的に応じた)ある程度の精度で分離しようと思えば、時刻tの状態からの変化が比較的小さいとみなせる時間的範囲(これを時刻tの近傍と呼ぶ)において前記複数の時刻を設定するのがよい。一例として、時刻t及び時刻tと隣接する前または後ろの時刻の透視画像を利用する態様をあげることができる。整合性の評価は、例えば、適当に定義された誤差の大きさ定量的に計算し、適当な値と比較することにより行なうことができる。
 そして、更新部では、前記計算部の評価結果に基づいて整合性を高めるように、言い換えれば分離の妥当性を高めるように、分布情報の更新を行なう。
 なお、分布情報は、透視画像内における対象部分あるいは背景部分の分布を表す情報である。分布情報は、例えば、空間的な位置情報や、各位置における輝度情報などによって表すことができる。透視画像をデジタル処理する場合、位置情報は、例えば、処理対象として着目する各領域(一画素からなってもよいし複数画素からなってもよい)における存在情報とすることができるし、輝度情報は、各領域を代表する輝度値の情報とすることができる。位置情報は、線形として扱われる場合には、例えば各領域の画素位置についての係数をもつ係数行列として表現することができる。もちろん位置情報を非線形として扱うことも可能であり、例えば、対象部分が拡大あるいは縮小を行なうような場合には、非線形の関数として表現することも有効と考えられる。
 対象部分と背景部分とは、統計的な解析を行なって分離されてもよい。例えば、対象部分及び背景部分の位置、形状、輝度のうち少なくとも一つの要素に対して、対象部分と背景部分との統計的独立性の有無・度合いなどを計算することにより、分離の妥当性を評価することができる。このような解析の例としては、独立成分分析を挙げることができる。
 (2)また、第2の案として、L(Lは自然数)個の対象部分と前記対象部分以外の背景部分との透過的な重ね合わせにより表現される透視画像を処理する信号処理装置または画像処理装置であって、ある時刻tにおいて、前記対象部分に関し透視画像を取得する透視画像取得部と、上記透視画像取得部で取得された上記透視画像のM×N(M,Nは自然数)個の領域に着目し、前記時刻tと前記時刻tよりも過去のK(Kは自然数)個の時刻t1,t2,…,tKとにおいて取得した(K+1)個の透視画像の(K+1)×M×N個の前記着目する領域が有する輝度値からなる(K+1)時刻分の観測ベクトルb(t,t1,t2,…,tK)を生成する観測ベクトル生成部と、前記着目する領域中の前記L個の対象部分それぞれに起因する輝度値からなる変数ベクトル(以下、対象変数ベクトルIa1,Ia2,…,IaLという)と前記着目する領域中の前記背景部分が有する輝度値からなる変数ベクトル(以下、背景変数ベクトルIという)とを、前記観測ベクトルbと前記時刻tよりも過去のある時刻tpにおける透視画像の前記着目する領域に関するL個の変数ベクトルIa1(tp),Ia2(tp),…,IaL(tp)とからなる対象変数ベクトルIa(tp)及び背景変数ベクトルIb(tp)を用いてそれぞれ更新する変数ベクトル更新部と、上記観測ベクトルbと、前記変数ベクトル更新部により更新された前記対象変数ベクトルI及び背景変数ベクトルIからなる変数ベクトルIと、前記変数ベクトルIを用いて推定された前記時刻tにおける前記L個の対象部分及び背景部分の位置と前記K個の過去の時刻t1,t2,…,tKにおいて同様に推定されたそれぞれの時刻における前記L個の対象部分及び背景部分の位置とに関する(K+1)MN×(L+1)MNの係数行列Aとから、fを決定論的、もしくは確率論的な関数として、式
Figure JPOXMLDOC01-appb-M000002
により定義される評価値PIの評価を高めるように前記変数ベクトルIを計算する変数ベクトル計算部と、をそなえ、上記透視画像取得部での前記透視画像の取得、上記観測ベクトル生成部での前記観測ベクトルbの生成、上記変数ベクトル更新部での前記変数ベクトルIの更新、及び、上記変数ベクトル計算部での前記変数ベクトルIの計算を、時刻をずらしながら行なうように構成された、信号処理装置または画像処理装置を用いることができる。
 (3)さらに、第3の案として、対象部分と前記対象部分とは異なる背景部分との透過的な重ね合わせにより表現される透視画像を処理するコンピュータに対し、複数の時刻について、前記対象部分を含む前記透視画像を取得する透視画像取得機能と、前記複数の時刻のうちのある時刻tにおける前記対象部分に起因する前記透視画像の成分及び前記背景部分に起因する前記透過画像の成分を、前記複数の時刻における前記対象部分及び前記背景部分の分布情報の推定値と関連づけて計算し、前記複数の時刻における前記透視画像との整合性を評価する計算機能と、前記計算機能による評価結果に基づいて、前記時刻tにおける前記対象部分又は前記背景部分の分布情報の推定値を更新する更新機能とを実現させるための信号処理プログラムを用いることができる。
 (4)また、第4の案として、上記プログラムを記録したコンピュータ読み取り可能な記録媒体を用いることができる。
 (5)さらに、第5の案として、上記の各処理を行なうステップを有する信号処理方法を用いることができる。
 上記本発明によれば、患部などの対象部分の位置,形状及び大きさなどを高精度に計測(推定)することが可能となる。
 また、金属性マーカを刺入することなく安全かつ低負担で、患部などの対象部分の位置,形状及び大きさなどを高精度に計測(推定)することが可能となる。
 さらに、放射線治療において、患部への正確な放射線照射を可能とすることが可能となる。
放射線治療装置の一例を示す図である。 図1に示す制御部の機能ブロック図である。 オプティカルフローの一例を示す図である。 ブロックマッチング(BM)法の一例を示す図である。 BM法により計測した患部位置と実際の患部位置との誤差を示す図である。 BM法の概略を示す図である。 テンプレートマッチング(TM)法の概略を示す図である。 透視画像と反射画像との相違点を説明する図である。 本例の信号処理方法の一例を示す図である。 本例の信号処理方法の一例を示す図である。 本例の信号処理方法の一例を示す図である。 本例の信号処理方法の一例を示す図である。 (A)は対象(腫瘍)部分の正解画像の一例を示す図であり、(B)は背景部分の正解画像の一例を示す図であり、(C)はファントムデータの透視画像の一例を示す図である。 ファントムデータの透視画像の一例を示す図である。 図14に示す透視画像から分離された対象部分を示す図である。 図14に示す透視画像から図15に示す対象部分が分離された背景部分を示す図である。 (A)~(C)は透視画像から分離された対象画像と背景画像とを比較した図である。 症例データの透視画像の一例を示す図である。 図18に示す透視画像から分離された対象部分を示す図である。 図18に示す透視画像から分離された背景部分を示す図である。 変形例の信号処理方法の一例を示す図である。 放射線治療装置の一例を示す図である。
 以下、図面を参照して本発明の実施形態を説明する。ただし、以下に示す実施形態は、あくまでも例示に過ぎず、以下に示す実施形態で明示しない種々の変形や技術の適用を排除する意図はない。即ち、本実施形態は、その趣旨を逸脱しない範囲で種々変形(各実施形態を組み合わせる等)して実施することができる。
 〔1〕一実施形態
 (1.1)放射線治療装置の構成例
 図1に放射線治療装置を例示する。この図1に示す放射線治療装置は、例示的に、制御部1と、放射線照射装置90と、計測用放射線照射部305と、センサ部306とをそなえ、例えば、ベッド230上に位置する患者の患部260に対して治療用の放射線を照射することができる。
 このとき、患部260の位置,形状及び大きさなどが、患者の呼吸運動などにより時間的に変動する場合がある。
 そこで、本例では、計測用放射線照射部305により透視画像撮影用の放射線(たとえば、X線など)を患部260周辺に照射し、患部260を通過した透視画像撮影用の放射線をセンサ部306により検知することにより、患部260周辺の放射線透視画像(以下、単に透視画像ともいう)を撮影し、制御部1へ入力する。そして、制御部1が、上記透視画像における患部260などの対象部分の位置,形状及び大きさなどの時間的な変動を計測(推定)する。また、制御部1は、当該計測結果に基づき、放射線照射装置90から照射される治療用の放射線の照射位置及び範囲を制御してもよい。このとき、放射線照射装置90から照射される治療用放射線を計測用放射線の代替えとして、患部260を通過した放射線によりElectronic Portal Imaging Device (EPID)などを用いて患部260の透視画像を撮影・計測(推定)してもよい。あるいは両者を融合して撮影・計測(推定)してもよい。さらに、制御部1はベッド230を制御して患部260と放射線の照射位置とを合わせてもよい。
 また、撮影方向の異なる透視画像を複数融合することにより、単一撮影面の2次元だけでなく3次元の位置、形状及び大きさとそれらの時間的変動(4次元)を同時に計測(推定)することが可能となる。このとき、複数の透視画像は同種の透視画像装置によって撮影されてもよいし、異種の装置によって撮影されてもよい。さらに、あらかじめ患部260の動きを撮影した4D-CTデータと透視画像データとを融合することによって、動きのある3次元計測を行なってもよい。
 なお、以下に説明する例では、基本的にX線シミュレータにより取得される透視画像を用いて患部260の位置,形状及び大きさなどの計測を行なった。なお、X線シミュレータでは、毎秒2フレームの透視画像を撮影することができるが、例えば、On Board Imager(OBI)では、毎秒15フレームの透視画像を撮影することができ、より高精度な計測を行なうことができる。
 (1.2)信号処理装置の一例
 ここで、制御部1は、図2に示すように、例えば、信号処理部2と、放射線照射制御部3とをそなえる。
 信号処理部(信号処理装置)2は、計測用放射線照射部305及びセンサ部306により撮影された透視画像に基づいて、時間的に変動する患部260の位置などを計測(推定)する。
 このため、信号処理部2は、例えば、透視画像取得部4と、観測ベクトル生成部5と、変数ベクトル更新部6と、変数ベクトル計算部7と、推定部8とをそなえる。なお、透視画像取得部4,観測ベクトル生成部5,変数ベクトル更新部6,変数ベクトル計算部7及び推定部8の各動作については、後記(1.6)にて詳述するが、ここでは各部の動作概要について説明する。
 透視画像取得部4は、計測用放射線照射部305及びセンサ部306により撮影された透視画像を取得する。前記透視画像は、ある時刻tを含む複数の時刻において取得され、透視画像には、位置,形状及び大きさなどの測定(推定)対象である対象部分が一般にL(Lは自然数)個含まれる。つまり、透視画像は、L個の対象部分と当該対象部分以外の背景部分(ただし、対象と背景との間に本質的な違いはないので、(L+1)個の対象と考えても良いし、いくつかの対象を背景と考えて複数の背景として扱っても良い)との透過的な重ね合わせにより表現されているといえる。
 また、観測ベクトル生成部5は、透視画像取得部4により取得された時刻tにおける透視画像のM×N(M,Nは自然数)個の画素(アナログ画像においてはデジタル化した画素)もしくは複数画素の何らかの代表値からなるある領域に着目して透視画像を記憶し、当該MN個の着目する領域が有する輝度値(以下、濃度ともいう)と、その時刻tよりも過去のK(Kは自然数)個の時刻t1,t2,…,tKにおいて前記と同様に取得・記憶したK個の透視画像の、あわせて(K+1)MN個の着目する領域からなる観測ベクトルb(t,t1,t2,…,tK)を生成する。なお、透視画像の記憶は上記透視画像取得部4で行なってもよい。
 変数ベクトル更新部(更新部)6は、観測ベクトル生成部5により生成された領域中のL個の対象部分それぞれに起因する輝度値からなる変数ベクトルIa1,Ia2,…,IaLと前記生成された領域中の背景部分が有する輝度値からなる変数ベクトルIとを、前記観測ベクトルbと前記時刻tよりも過去のある時刻tpにおける透視画像の前記着目する領域に関するL個の変数ベクトルIa1(tp),Ia2(tp),…,IaL(tp)とからなる対象変数ベクトルIa(tp)及び背景変数ベクトルIb(tp)を用いてそれぞれ更新する。
 また、変数ベクトル計算部(計算部)7は、前記観測ベクトル生成部5により生成された上記観測ベクトルbと、前記変数ベクトル更新部6により更新された前記対象変数ベクトルI及び背景変数ベクトルIからなる変数ベクトルIと、前記変数ベクトルIを用いて推定された前記時刻tにおける前記L個の対象部分及び背景部分の位置と前記K個の過去の時刻t1,t2,…,tKにおいて同様に推定されたそれぞれの時刻における前記L個の対象部分及び背景部分の位置とに関する(K+1)MN×(L+1)MNの係数行列(関数)Aとから、fを決定論的、もしくは確率論的な関数として、式
Figure JPOXMLDOC01-appb-M000003
により定義される評価値PIを、たとえば最適にするなどの適切な値にするように前記変数ベクトルIを計算する。
 推定部8は、前記変数ベクトル更新部6によりある初期値が設定された対象変数ベクトルIと前記変数ベクトル更新部6により更新された対象変数ベクトルIとに基づいて、前記対象部分の時間的な位置変動を推定してもよいし、あるいは、前記変数ベクトル更新部6により更新される前の対象変数ベクトルIと前記変数ベクトル更新部6により更新された後の対象変数ベクトルIとに基づいて、前記対象部分の時間的な位置変動を推定してもよい。また、更新された対象変数ベクトルIaに基づいて、対象の形状及び大きさなどを推定してもよい。
 そして、上記の信号処理部2は、上記透視画像取得部4での透視画像の取得、上記観測ベクトル生成部5での前記観測ベクトルbの生成、上記変数ベクトル更新部6での前記変数ベクトルIの更新、及び、上記変数ベクトル計算部7での前記変数ベクトルIの計算を、時刻をずらしながら行なうように構成される。
 信号処理部2により計測された患部260などの対象部分の位置,形状及び大きさなどの情報は、放射線照射制御部3へ出力される。
 放射線照射制御部3は、信号処理部2からの前記情報に基づいて、放射線照射装置90から照射される治療用の放射線の照射位置及び照射範囲、あるいはベッド230の位置などを制御する。
 (1.3)一般的な位置推定手法の検証
 まず、上記のような透視画像について、図3に示すように、一般的な位置推定方法として知られるオプティカルフロー抽出を試みた。これは、異なる時間において撮影された透視画像の各フレーム(動画像フレームともいう)間における対象部分(例えば、患部260など)Tの各画素について移動量ベクトル(u,v)を抽出するものである。なお、以下では簡単のため対象部分の個数を1(L=1)として説明する。
 図3は、時刻t(>0)における対象部分Tと、離散的な次の時刻(t+1)における対象部分Tとについて、上記移動量ベクトル(u,v)を求める例である。ここで、次の時刻(t+1)はtよりも未来の任意の時刻でも良いが、以下では簡単のため単に(t+1)とする。
 まず、オプティカルフロー抽出法の一例として、一般的なブロックマッチング(BM)法を用いる例について説明する。なお、BM法をはじめ、勾配法などの一般的なオプティカルフロー抽出法では、対象部分の各画素の対応する輝度値が、対象部分の移動によっても変化しない(一定である)ことを前提条件としている。
 BM法は、ある時刻tで撮影されたフレーム(前フレームともいう)において対象部分のある注目点を含むある大きさの領域をテンプレートのブロックとして、時刻(t+1)で撮影された次のフレーム(次フレームともいう)のある決められた範囲(画像すべてでも良い)を全探索することにより、前フレームのテンプレートブロックとの間で差分などの評価関数を最適化する次フレームのブロックの注目点とブロック内の相対位置が同じ点を対応点とすることで、フレーム間での対象部分を構成する各点の移動変位を抽出する手法である。
 具体的には、例えば、図4に示すように、まず、時刻tに撮影されたフレーム(前フレーム)において、対象部分を含むM×N(図4に示す例では、M=N=3)のあるテンプレートブロックに対して、時刻(t+1)に撮影されたフレーム(次フレーム)における各ブロックのうち、前フレームにおけるあるブロックと差分評価関数等を最小とするブロックを探索し、当該ブロック間で移動量ベクトル(u,v)を算出する。
 上記算出は、例えば、以下のような計算により行なわれる。
 例えば、時刻tに撮影されたフレームにおける対象部分の位置(m,n)におけるある画素の輝度値をI(m,n,t)とし、時刻(t+1)に撮影されたフレームにおける抽出対象の位置(m,n)における前記画素の輝度値をI(m,n,t+1)とした場合、前フレームにおけるある画素の輝度値I(m,n,t)と、次フレームにおいて当該画素の位置から移動量ベクトル(u,v)だけ移動した画素の輝度値I(m+u,n+v,t+1)との最小二乗誤差Mean Square Error(MSE)は次式(a)で定義される。
Figure JPOXMLDOC01-appb-M000004
 ここで{m,n}は、上記M×Nのブロック内に存在するすべての画素位置の集合である。BM法では、たとえばこのMSEを最小にするフロー(U,V)が、次式(b)により、各ブロックに対してそれぞれ演算される。
Figure JPOXMLDOC01-appb-M000005
 実際にBM法を用いて患部260の位置を計測(推定)したところ、患部260の実際の位置と推定した位置との差(推定誤差)は、図5に示すように、約1.2[mm]~約2.4[mm]となった。ただしこの方法は、対象部分の輪郭が明瞭な特定の一部分を用いるなど、高コントラスト部分を用いることで精度の向上が期待できるため、対象の患部260の腫瘍部分のうち、目視によりコントラストが高い輪郭部分を注目点として手動設定した。また、ブロックサイズは正方のN×Nとした。なお、図5中、横軸はブロックサイズ(N)の大きさを表す。
 このように、患部260の位置推定にBM法を用いた場合、上記のような精度向上の工夫をしたにも拘らず、臨床上要求される精度(例えば、推定誤差が1[mm]以内)を達成できない。
 (1.4)一般的な位置推定手法の問題点
 上記のように、BM法では、所望の推定精度(誤差が1[mm]以内)を達成できないが、この原因としては、主に次の2つが考えられる。
 第1に、透視画像は非常に不明瞭であり、同じ対象であってもノイズなどによりフレーム間で輝度値が異なる可能性があるということ。
 第2に、透視画像において、患部260の輝度値は、患部260などの対象部分そのものの輝度値と背景部分の輝度値との積算値により表されるため、上記オプティカルフロー抽出法の前提条件(即ち、対象部分の輝度値が対象部分の移動によっても変化しない)が成り立たない。この条件は、位相限定相関法などのオプティカルフロー抽出以外の他の一般的な方法でも前提としているため、当該条件の不成立は、精度低下を招く、より本質的な原因となっている。
 そこで、上記の各問題を解決する方法について以下に述べる。
 (1.5)テンプレートマッチング
 BM法では、図6に例示するように、前フレームと次フレームとの間で対象部分の差分などの評価値を算出するため、前フレーム及び次フレームのうち少なくとも一方の透視画像が不明瞭であったり、同一対象部分であってもノイズなどによりフレーム間で輝度値が異なる場合、評価値の算出精度が低下し、位置推定精度が低下する。
 そこで、BM法ではなく、例えば、テンプレートマッチング(TM:Template Matching)法を用いることができる。
 TM法は、図7に例示するように、対象部分のテンプレートを予め取得しておき、当該テンプレートと最も相関の高い範囲を、各時刻におけるフレームにおいてそれぞれ探索する方法である。なお、TM法に用いるテンプレートは、例えば、ある時刻におけるフレーム中、所定値以上の輝度値を有する範囲をテンプレートとして定めることで取得することができる。または、他の計測装置(CTやMRI(magnetic resonance imaging)など)により計測された患部260の位置,形状及び大きさなどに基づいて、テンプレートを定めてもよい。
 TM法では、たとえばテンプレートが有する輝度値Itemp(m,n)と、時刻tのフレームにおける対象部分の輝度値I(m+u,n+v,t)との最小二乗誤差(MSE)を次式(c)により定義する。
Figure JPOXMLDOC01-appb-M000006
 そして、TM法では、たとえば上記の式(c)におけるMSEを最小にするフロー(U,V)を、次式(d)により算出する。
Figure JPOXMLDOC01-appb-M000007
 このように、TM法では、ノイズの影響を受けない対象部分のテンプレートと各フレームにおける対象部分との評価値に基づき、対象部分の位置を計測(推定)するので、その推定結果が、比較対象両者がノイズの影響を受ける可能性のあるBM法に比して、透視画像の不明瞭さに影響されにくい。
 また、BM法では、対象部分の細部にわたる濃度分布が評価値に影響するが、TM法においては、例えば、濃度分布一定のテンプレートを用いることにより、対象部分の細部の相違にはあまり影響されない効果が期待できる。
 この他、ノイズの影響を低減するため、式(c)にノイズ項を付加して統計的な処理を行なうことも有効である。
 (1.6)濃度分布の分離
 また、図8(A)に示すように、透視画像においては、患部260などの対象部分の見かけ上の輝度値は、対象部分本来の輝度値と、患部260の周辺臓器などの背景部分の輝度値との積算値で定まる。
 そのため、たとえ対象部分の一部分に着目しても、当該一部分に重なる背景部分の輝度値が時間的に異なっていれば、当該一部分の見かけ上の輝度値は時間的に一定とならない。従って、一般的には、透視画像において、前述した従来のオプティカルフロー法の前提条件が成立しない。
 一方、図8(B)に示すように、透視画像でない通常の画像(反射画像という)において、対象部分の輝度値は、背景部分の輝度値に応じて変化しない。つまり、反射画像において対象部分の一部分に着目した場合、当該一部分に重なる背景部分の輝度値が時間的に異なっていても、当該一部分の見かけ上の輝度値は一定である。従って、一般的に、反射画像においては、前述した従来のオプティカルフロー法の前提条件が成立する。
 そこで、本例では、図9に例示するように、透視画像における見かけ上の(観測された)輝度値I(m,n,t)を、時間的に移動する対象部分本来の輝度値I(m,n,t)と、対象とは別の動きをする背景部分の輝度値I(m,n,t)とに分離する。なお、以下では引き続き簡単のため対象部分の個数を1(L=1)とし、背景の移動変位を0(つまり、静止している)として、単にI(m,n)と記す。
 透視画像から分離された対象部分の輝度値に関しては、オプティカルフロー法の前提条件が成立するので(本例では、この前提条件が成り立つような分離を行なう)、透視画像において観測された輝度値I(m,n,t)ではなく、時間的に移動する対象部分本来の輝度値I(m,n,t)について、TM法などにより位置推定を行なうことにより、反射画像と同様に対象部分の位置を高精度に計測(推定)することが可能となる。また、推定された位置情報を用いて、Iaの分離精度を高める事が可能なため、各処理を繰り返すことにより対象部分の位置だけでなく、形状及び大きさなどを高精度に計測(推定)することが可能である。
 以下に、本例の分離方法について、図10~図12を用いて説明する。
 まず、図10に例示するように、例えば、時刻tにおける透視画像のある領域(注目領域(ROI:Region of Interest)ともいう)の輝度値ベクトル(以下、観測ベクトルともいう)I(m,n,t)は、対象部分の輝度値ベクトルI(m,n,t)と、背景部分の輝度値ベクトルI(m,n)とにより、次式(e)のように表すことができる。
Figure JPOXMLDOC01-appb-M000008
 なお、上記の式(e)において、Iamn,(m,n)={(2,2),(2,3),(3,2),(3,3)}は、透視画像のROI中の各ブロックの位置座標(m,n)における対象部分本来の輝度値を示した変数ベクトルである。また、Ibmn,(m,n)={(2,2),(2,3),(3,2),(3,3)}は、透視画像のROI中の各ブロックの位置座標(m,n)における背景部分の輝度値を示した変数ベクトルである。
 この図10から分かるように、上記式(e)で表される方程式の解は、
Figure JPOXMLDOC01-appb-I000009
であるが、式(e)で表される方程式は未知変数の数に対して独立な方程式の数が不足している不定であり、対象部分の変数ベクトルI(Ia22~Ia33)及び背景部分の変数ベクトルIの分布(Ib22~Ib33)は一意に定まらない。これは、1枚の透視画像からその濃度を構成している要素部分を一般には分離できないことを意味している。
 そこで、図11に例示するように、時刻tにおいて取得された観測フレームと別の時刻(簡単のため引き続き(t+1)とする)において取得された観測フレームを時間的に統合して、各ROIの観測ベクトルIを、対象部分の変数ベクトルIと、背景部分の変数ベクトルIとにより、次式(f)のように表す。
Figure JPOXMLDOC01-appb-M000010
 上記の式(f)において、左辺左側の行列は係数行列(関数)Aであり、例えば、異なる時刻t,(t+1)における対象部分の位置により決定される行列である。また、左辺右側の行列は、時刻tにおけるROIの対象部分及び背景部分の濃度を示す変数ベクトルIである。また、右辺の行列は、異なる時刻t,(t+1)において取得された透視画像のROIにおいて観測される輝度値を示す観測ベクトルbである。なお、本例では、透視画像を4×4の16画素もしくは複数画素の何らかの代表値16個から構成されるものとし、ROIを2×2の4(M=2, N=2よりMN=4)画素で規定しているので、2時刻分をまとめた係数行列は、K=1, L=1より(K+1)MN=8, (L+1)MN=8だから8×8の行列となり、変数ベクトルIは8×1の行列となっているに過ぎず、透視画像の画素(もしくは代表値)数やROIの画素(もしくは代表値)数、ならびに対象部分の個数、統合する過去の時刻数に応じて各行列のサイズはそれぞれ変化する。
 ここで、係数行列Aがフルランクであれば、当該係数行列Aで関連付けられる変数ベクトルIは一意に求まる。
 しかし、実際には、異なる時刻t、時刻(t+1)間において、対象部分がどのように移動したのかは不明である、即ち、係数行列A(特に、係数行列Aの5,6,7,8行目のそれぞれ1,2,3,4列目を抽出した4×4の部分行列{aij}, i=5,6,7,8, j=1,2,3,4, ただし、aijは行列Aのi行j列目の要素である。なお、背景は静止しているものとしているが、背景も移動する場合はj=1,…,8である。)が不明であるため、本例では、係数行列Aを推定する。
 また、係数行列Aで関連付けられる変数ベクトルIの各要素は、例えば、既述のテンプレート作成方法により作成したテンプレートが有する輝度値を初期値として設定することができるが、あくまで仮定した値であるから一般に実際の値とは異なる。
 そこで、本例では、上記の係数行列A及び変数ベクトルIを、ある初期値から開始して相互の推定値を用いた更新を繰り返すことで推定精度を高める。
 図12は本例の動作について示したフローチャートの一例である。
 この図12に例示するように、まず、信号処理部2は、初期設定として、時刻パラメータtに「0」を代入し、透視画像における注目領域(ROI)を設定する。そして、変数ベクトル更新部6が、変数ベクトルIの初期値として、例えば、対象部分のテンプレートの輝度値ベクトルI templateを設定する(ステップA1)。
 ここで、ROIは、患部260などの対象部分を含むように設定されるのが望ましい。例えば、対象部分を別に撮影したCTやMRI画像などと対象部分の移動範囲の概略とを元に、常に対象部分が含まれるような領域を自動で、あるいは目視による手動で、またはあらかじめ取得したある時刻の透視画像において、所定値以上の輝度値を有する部分(範囲)を対象部分であると仮定することにより、当該部分が動く範囲を見積もってそれを含む領域をROIとして設定することができる。
 また、テンプレートの輝度値ベクトルI templateについては、変数ベクトル更新部6が、例えば、対象部分を別に撮影したCTやMRI画像などを元に自動で、あるいは目視による手動で対象部分の輪郭を大まかに抽出、またはある時刻に取得された透視画像において、所定値以上の輝度値を有する部分(範囲)を対象部分であると仮定することにより、当該部分を対象部分のテンプレートとして設定することができる。さらに、抽出部分の輝度値を適当な一定値とすることで、あるいは当該テンプレートが有する明るさに関する情報に基づいて、I templateを算出することができる。
 次に、透視画像取得部4が、時刻tにおいて、対象部分に関し透視画像を取得する。即ち、透視画像取得部4及び観測ベクトル生成部5により、時刻tにおいて取得された透視画像について、ステップA1で設定されたROIにおける輝度値ベクトルI(t)が取得される(ステップA2)。
 そして、変数ベクトル更新部6が、次のステップA3~A8に基づいて、対象部分の移動量ベクトル(u,v)を推定するために用いられる、変数ベクトルIの推定値である
Figure JPOXMLDOC01-appb-I000011
の更新処理を行なう(図12中の破線部分を参照)。なお、図12中、各変数に付されたハット記号は、当該変数が推定値であることを示すものである。また、行列tAは行列Aの転置行列を表す。なお、以下の説明では、表示上の制約から、ハット記号の代わりにダッシュ記号「′」を用いることがある。
 当該更新処理においては、まず、変数ベクトル更新部6により、時刻パラメータtが、t=0を満たすかどうかが判定される(ステップA3)。
 ここで、t=0を満たすと判断された場合(ステップA3のYesルート)、変数ベクトル更新部6は、変数ベクトルI′の初期化を行なう(ステップA4)。ステップA4でのI′の初期化は、例えば、観測ベクトル生成部5により生成された、時刻t=0における観測ベクトルI(0)から、上記テンプレートの輝度値ベクトルI templateを減算し、当該減算結果がI′(0)に設定されることで行なってもよい。
 また、観測ベクトル生成部5は、時刻t=0における観測ベクトルI(0)をb(0)とする(ステップA5)。
 一方、t=0を満たさないと判断された場合(ステップA3のNoルート)、観測ベクトル生成部5は、例えば、上記取得された透視画像のROIをM×N個の画素(もしくは代表値)領域とし、当該MN個のROIが有する輝度値と、時刻(t-1)において前記と同様に取得済の透視画像のROIのMN個の画素(もしくは代表値)の輝度値に基づいてあわせて2MN個の画素(もしくは代表値)の観測ベクトルb(t,t-1)を生成する(ステップA6)。なお、輝度値については、例えば、観測ベクトル生成部5が、透視画像のROIにおいて、各ブロックでの単位面積あたりの明るさを測定し、当該測定結果に基づいて算出するようにしてもよい。そして、変数ベクトル更新部6は、I′の更新を行なう(ステップA7)。ステップA7でのI′の更新は、例えば、変数ベクトル計算部7により前の時刻で計算された変数ベクトルI′の値が、次の時刻における変数ベクトルI′に設定されることで行なわれる。
 次に、変数ベクトル更新部6は、変数ベクトルI′の更新を行なう(ステップA8)。ステップA8でのI′の更新は、例えば、観測ベクトル生成部5により生成された、時刻tにおける観測ベクトルI(t)から、ステップA4またはA7において初期化または更新された変数ベクトルI′を減算し、当該減算結果がI′に設定されることで行なわれる。
 そして、推定部8が、例えば、ステップA1において設定されたI templateとステップA8において更新された変数ベクトルI′との間でTM法や勾配法を適用することにより、対象部分の移動量ベクトルの推定値(u´(t),v´(t))を算出(推定)する(ステップA9)。このとき、対象部分を、対象部分を構成するすべての点が同じ移動量をもつ剛体と仮定することにより、推定処理を高速化してもよいし、すべての点もしくは一部の点が他の点とは異なる移動量をもつ非剛体として厳密に推定してもよい。以下では簡単のため、対象を剛体と仮定して説明する。なお、推定部8は、例えば、ステップA8において更新される前の変数ベクトルI′とステップA8において更新された後の変数ベクトルI′との間でBM法や勾配法を適用することにより、対象部分の移動量ベクトルの推定値(u´(t),v´(t))を算出(推定)してもよいし、Ia templateもしくはステップA8において更新された後の変数ベクトルIa’との位相限定相関法やパーティクルフィルタを用いた方法などで算出(推定)してもよい。
 次に、変数ベクトル計算部7が、推定部8により推定された移動量ベクトルの推定値(u´(t),v´(t))に基づいて、係数行列Aを推定・更新する。即ち、変数ベクトル計算部7(関数更新部)は、初期値が設定された変数ベクトルIと更新された変数ベクトルIとに基づいて、係数行列(関数)Aを推定・更新する機能を有する。あるいは、変数ベクトル計算部7は、更新される前の変数ベクトルIと前記更新された後の変数ベクトルIとに基づいて、係数行列Aを推定・更新する機能を有する。
 また、t=0でない場合は変数ベクトル計算部7が、たとえば次式(g)により定義される誤差eの各要素の2乗を最小にするような、I(t)を計算により求める(ステップA10)。
Figure JPOXMLDOC01-appb-M000012
 ここで、行列Aの真値が求まり、かつフルランクであれば、誤差e=0を満足する連立方程式の厳密解Iを行列Aの逆行列を計算することで容易に求める事ができる。しかし、上述のように実際には行列Aは推定値であり、仮にフルランクになった場合の厳密解が求まったとしても、それが分離したいIa及びIbに一致する保証はない。そこで、本例では厳密に解を求めるのではなく、最適化手法などを用いて分離濃度を推定する。
 誤差eの2乗の最小化方法には、例えば、最急降下法や共役勾配法、Newton法などの最適化法をはじめ、人工神経回路網やファジィ理論、進化的計算理論などのいわゆる計算知能的手法も用いることができる。また、誤差関数に統計的なノイズ項を付加して、確率的な最適化法を用いることもできる。これはノイズの多い不明瞭な画像では有効な方法である。さらに、本例で対象とする腫瘍部分に起因する画像は、対象だけが存在し対象部分以外の輝度が0であるものが多いため、推定されたI’aに対し、輝度値がある値よりも小さな画素の値を閾値処理などによりノイズ除去してもよい。これにより特に最適化の過程での推定誤差の蓄積軽減効果が期待できる。
 例えば、最急降下法を用いることにより誤差eの各要素の2乗を最小にするようなI(t)を求める場合、以下の式(h)が用いられる。
Figure JPOXMLDOC01-appb-M000013
 上記の式(h)において、b(t,t-1)=t[tI(t-1),tI(t)]=b(I)であり、ηは変更割合を表すある正の係数である。
 即ち、変数ベクトル計算部7は、前記推定した係数行列Aに基づいて最適化法を用いることにより、誤差eのような何らかの評価指標を適切にするような、変数ベクトルIを計算するのである。
 そして、信号処理部2は、時刻tの値をインクリメントし(ステップA11)、ステップA2~ステップA11の処理を繰り返し実行する。
 即ち、本例の信号処理部2は、透視画像取得部4での前記透視画像の取得、観測ベクトル生成部5での前記観測ベクトルbの生成、変数ベクトル更新部6での前記変数ベクトルI及びIの更新、及び、変数ベクトル計算部7での前記変数ベクトルI及びIの計算を、時刻をずらしながら行なうように構成される。
 従って、変数ベクトル更新部6は、変数ベクトル計算部7により計算された変数ベクトルI及びIを用いて、変数ベクトルI及びIをそれぞれ更新するのである。
 上述したように、本例の信号処理方法によれば、透視画像から対象部分と背景部分とを分離することができ、分離された対象部分について位置推定を行なうことにより、対象部分の位置を高精度に測定(推定)することが可能となる。また、分離により対象部分の形状や大きさなどを高精度に測定(推定)することも可能となる。
 (1.7)実験例
 はじめに正解の分離濃度が分かっているファントムデータを用いた実験結果について説明する。ファントムデータとは、ROIのサイズが100×100の対象腫瘍がない図13(B)のような透視画像を背景とし、当該背景に図13(A)のような人工的に移動する腫瘍画像を積算した図13(C)に示すような透視画像データである。図14に時間発展に伴う上記ファントムデータのフレームの変化例を示す。図14に示す例では、合計110フレームの中から適当に選択した12フレームの透視画像が、時系列順に、左上から右下に向かって配置されている。腫瘍はROIの中央付近を主に縦方向に移動しており、実際の肺腫瘍の動きを模擬したものとなっている。
 この図14に示す各フレームについて、本例の信号処理方法を適用したところ、図15及び図16に示すように、正解とはかなり異なる初期値を手動で適当に設定した場合でも、時間の進展ともに対象部分と背景部分とを精度よく分離することができた。特に、図15に示すように、目視で適当に手動設定した初期値から、時間発展に伴いほぼ図13(A)に示す腫瘍画像を分離抽出できていることが分かる。初期値の誤差がROIサイズ10000個のうち300画素以上で、濃度差も1画素あたり20だったのに対し,最終フレームでは誤差は0画素で、1画素あたりの濃度差も平均約0.1以下まで減少させることに成功しており、高精度に分離できていることが示された。対象腫瘍の形状そのものは完全に分離できていることから、残りの誤差も対象腫瘍の形状内の濃度値の僅かな誤差のみである。この事は、図16に示すように腫瘍を分離した背景画像が図13(B)に一致していることからも確認できる。初期値設定と更新とを繰り返して得られた分離結果を正解画像と比較した例を図17に示す。初期画像(B)では腫瘍部分の形状や大きさなどが正解画像(A)とは大きく異なり、背景画像にもその影響がはっきり見て取れるが、分離結果画像(C)では腫瘍の形状や大きさも正確に分離できている。また、濃度分布に多少の誤差が残っているが、本例の信号処理方法をさらに繰り返し適用することで、さらなる高精度な推定が可能であり、この程度の誤差であれば、背景画像にその影響を目視するのは困難である。
 また、このとき最終的に分離されたIaについて、位置推定手法の一例としてTM法を用いて対象部分の位置計測を行なったところ、その推定誤差はほぼ0となり正確な位置推定結果が得られた。同じファントムデータに従来のBM法を適用したところ、背景画像の高コントラストな肋骨輪郭などに影響され、平均誤差は約3[mm]となってしまったことからも、本信号処理方法の有効性が確認できる。
 次に、実際の症例データを用いた実験結果について説明する。図18に、ROIサイズ60×60の対象部分と背景部分とを分離する前の透視画像の一例を示す。図18に示す例では、第1フレームから第12フレームまでの連続した透視画像が、時系列順に、左上から右下に向かって配置されている。
 この図18に示す各フレームについて、本例の信号処理方法を適用したところ、図19及び図20に示すように、対象部分と背景部分とを分離することができた。
 ここで、図19に示す分離された対象部分について、位置推定手法の一例としてTM法を用いて対象部分の位置計測を行なったところ、その推定誤差は、0.61±0.18[mm]となった。また、対象部分の形状及び大きさについても高精度に測定することができた。ファントムデータが100フレーム分の更新だったのに対し、用いた症例データは12フレーム分しか使用できない制約があったため、ファントムデータに比べると精度はやや落ちるものの、臨床上十分な高精度計測が可能となっている。また、実際の臨床では治療照射前に1,2分程度の観測を行なうことで、100フレーム程度のデータを取得することが可能なため、症例データに対してもさらに高精度な分離と位置計測が可能である。
 これに対し、対象部分と背景部分とを分離することなく、図18に示す透視画像について、位置推定手法の一例としてTM法を用いて対象部分の位置計測を行なったところ、その推定誤差は、0.81±0.04[mm]であった。BM法よりも精度が向上しており、TM法の効果が確認できる。ただし、この場合のTM法は、BM法と同様に対象部分の輪郭が明瞭であるなど、高コントラストであれば高精度が期待できるため、目視によりコントラストが高い輪郭部分を注目点として試行錯誤的に手動設定し、最も推定誤差が小さくなるROIサイズを採用した。目視による手動設定は細心の注意を要する煩雑な作業であるため、多忙な臨床現場で行なうことは設定者の負担も大きい。したがって、対象と背景とを分離することなく高精度な計測を行なうことは難しい。
 以上から、本例の信号処理方法のように、透視画像から対象部分と背景部分とを分離し、分離された対象部分について位置推定を行なうことにより、異なる時刻において取得された透視画像から対象部分の位置,形状及び大きさなどを高精度に計測(推定)することが可能となる。
 また、本例の信号処理方法によれば、観測対象に金属性マーカを刺入することなく安全かつ低負担で、対象部分の位置,形状及び大きさなどを高精度に計測(推定)することができる。
 (1.8)変形例
 本例では、図21の丸付き数字1に例示するように、対象部分の少なくとも一部が、透視画像に設定されたROIの外へ移動した場合の信号処理方法の一例について説明する。
 例えば、図11に示す例では、時刻tにおいて、(m,n)=(2,2)及び(3,2)の位置にあった対象部分の一部が、時刻(t+1)において、(m,n)=(2,1)及び(3,1)の位置へ移動している。
 このような場合、ROIの他の部分(図11に示す例では、時刻(t+1)の透視画像において、(m,n)=(2,3)及び(3,3)の部分)に対応する変数ベクトルIの各要素の係数は0となる。(ただし、上記の式(f)を連立方程式として定式化するため、I′を構成する変数は一定とする。)
 上記の例では、時刻tにおける値を上付き表示でIij(t)=Iij tのようにも表すと
Figure JPOXMLDOC01-appb-I000014
となるが、一方で、
Figure JPOXMLDOC01-appb-I000015
となるため、矛盾が生じる。
 本例では、例えば、上記の矛盾を、
Figure JPOXMLDOC01-appb-I000016
とすることで回避することができる(図21の丸付き数字2参照)。
 しかしながら、このまま行列計算
Figure JPOXMLDOC01-appb-I000017
により解くと、誤った拘束条件
Figure JPOXMLDOC01-appb-I000018
により、
Figure JPOXMLDOC01-appb-I000019
が更新されない。
 その結果、当該拘束条件の下で、他の変数のみを計算することになり、
Figure JPOXMLDOC01-appb-I000020
が真値でなければ、他の変数は誤った値となる。
 そこで、誤った拘束条件の影響を排除するため、前述の最急降下法などの最適化手法を用いることができる。
 例えば、誤差eを次式(i)で定義し、当該誤差eについて最急降下法などの最適化手法を適用することにより、I´を計算することができる(図21の丸付き数字3参照)。
Figure JPOXMLDOC01-appb-M000021
 これにより、対象部分の一部がROIの外へ移動した場合であっても、より高精度に対象部分の濃度分布を算出することができる。結果、このような場合であっても、対象部分の位置,形状及び大きさなどをより高精度に計測(推定)することが可能となる。
 (1.9)その他
 また、上述した信号処理部2により推定された対象部分の位置に関する情報を用いて、前記対象部分に放射線を照射する放射線治療装置を実現することもできる。
 例えば、図22に例示する放射線治療装置において、患者の患部260を前記対象部分として、放射線を発生する放射線発生部(パルスモジュレータ250,クライストロン100,冷却装置110,導波管120,真空ポンプ130,電子銃140,加速管150,偏向電磁石160及びターゲット170)と、前記放射線発生部から発生した前記放射線の照射範囲を所望の形状に形成するコリメータ部(MLC210)と、前記患部の時間的な位置変動を推定する信号処理装置(信号処理部2)と、前記信号処理装置により推定された前記患部の位置に関する情報を用いて、前記放射線の照射位置と照射範囲とを算出し、前記算出結果を基に前記コリメータ部を駆動制御する駆動制御部(制御部1)と、をそなえた、放射線治療装置が考えられる。
 これにより、放射線治療において、患部への正確かつ連続的な放射線照射が可能となる。 
 なお、上述した信号処理部2としての機能は、コンピュータ(CPU,情報処理装置,各種端末を含む)が所定のアプリケーションプログラム(信号処理プログラム)を実行することによって実現されてもよい。
 即ち、上記プログラムは、対象部分と前記対象部分とは異なる背景部分との透過的な重ね合わせにより表現される透視画像中の前記対象部分の時間的な位置変動を推定する信号処理装置において、前記推定機能をコンピュータに実現させるための信号処理プログラムであって、前記信号処理プログラムが、ある時刻tにおいて、前記対象部分に関し透視画像を取得する透視画像取得機能と、上記透視画像取得機能により取得された上記の透視画像のM×N個の画素(アナログ画像においてはデジタル化した画素)もしくは複数画素の何らかの代表値からなるある領域に着目し、前記ある時刻tよりも過去のK(Kは自然数)個の時刻t1,t2,…,tKにおいて前記と同様に取得したK個の透視画像の、あわせて当該(K+1)MN個の前記着目する領域が有する輝度値からなる(K+1)時刻分の観測ベクトルb(t,t1,t2,…,tK)を生成する観測ベクトル生成機能と、前記着目する領域中の前記L個の対象部分それぞれに起因する輝度値からなる変数ベクトルIa1,Ia2,…,IaLと前記着目する領域中の前記背景部分が有する輝度値からなる変数ベクトルIとを、前記観測ベクトルbと前記ある時刻tよりも過去のある時刻tpにおける透視画像の前記着目する領域に関する変数ベクトルIa1(tp), Ia2(tp),…, IaL(tp) からなる対象変数ベクトルIa(tp)及び背景変数ベクトルIb(tp)とを用いてそれぞれ更新する変数ベクトル更新機能と、上記観測ベクトルbと、前記変数ベクトル更新機能により更新された前記変数ベクトルI及びIからなる変数ベクトルIと、前記変数ベクトルIを用いて推定された前記時刻tにおける前記L個の対象部分及び背景部分の位置と前記K個の過去の時刻t1,t2,…,tKにおいて同様に推定されたそれぞれの時刻における前記L個の対象部分及び背景部分の位置とに関する(K+1)MN×(L+1)MNの係数行列Aとから、fを決定論的、もしくは確率論的な関数として、式
Figure JPOXMLDOC01-appb-M000022
により定義される評価値PIを、たとえば最適にするなどの適切な値にするように前記変数ベクトルIを計算する変数ベクトル計算機能と、をそなえ、上記透視画像取得機能による前記透視画像の取得、上記観測ベクトル生成機能による前記観測ベクトルbの生成、上記変数ベクトル更新機能による前記変数ベクトルIの更新、及び、上記変数ベクトル計算機能による前記変数ベクトルIの計算を、時刻をずらしながら行なうように前記コンピュータを動作させる、信号処理プログラムの一例である。
 また、上記プログラムは、例えば、フレキシブルディスク,CD(CD-ROM,CD-R,CD-RWなど),DVD(DVD-ROM,DVD-RAM,DVD-R,DVD-RW,DVD+R,DVD+RWなど)等の一時的でない(non-transitory)コンピュータ読取可能な記録媒体に記録された形態で提供されうる。この場合、コンピュータはその記録媒体から信号処理プログラムを読み取って内部記憶装置または外部記憶装置に転送し格納して用いることができる。また、そのプログラムを、例えば、磁気ディスク,光ディスク,光磁気ディスク等の記憶装置(記録媒体)に記録しておき、その記憶装置から通信回線を介してコンピュータに提供するようにしてもよい。
 ここで、コンピュータとは、ハードウェアとOS(オペレーティングシステム)とを含む概念であり、OSの制御の下で動作するハードウェアを意味している。また、OSが不要でアプリケーションプログラム単独でハードウェアを動作させるような場合には、そのハードウェア自体がコンピュータに相当する。ハードウェアは、少なくとも、CPU等のマイクロプロセッサと、記録媒体に記録されたコンピュータプログラムを読み取るための手段とをそなえている。
 上記信号処理プログラムとしてのアプリケーションプログラムは、上述のようなコンピュータに、信号処理部2としての機能を実現させるプログラムコードを含んでいる。また、その機能の一部は、アプリケーションプログラムではなくOSによって実現されてもよい。
 なお、本実施形態としての記録媒体としては、上述したフレキシブルディスク,CD,DVD,磁気ディスク,光ディスク,光磁気ディスクのほか、ICカード,ROMカートリッジ,磁気テープ,パンチカード,コンピュータの内部記憶装置(RAMやROMなどのメモリ),外部記憶装置等や、バーコードなどの符号が印刷された印刷物等の、コンピュータ読み取り可能な種々の媒体を利用することもできる。
 また、上述した放射線治療装置及び信号処理部2の各構成及び各処理は、必要に応じて取捨選択してもよいし、適宜組み合わせてもよい。
 例えば、変形例では、対象部分がROIの外へ移動した場合の位置推定方法について説明したが、ROIの範囲を十分大きく設定して実際には推定しないダミー変数を用いるなどすることにより、当該事態を回避してもよい。ただし、ダミー変数を用いる場合、逆行列計算で解を計算することはできず、変形例で示した誤った拘束条件を修正しながら最適化手法を用いる方法などを用いる。
 また、本例の信号処理方法は、例えば、血管造影法により撮影される血管についての位置,形状及び大きさなどの推定や、透過性のある複数の微生物などが一部もしくは全部重なりあって運動している場合の顕微鏡映像における複数対象物のそれぞれの位置変動、形状及び大きさ変動などの推定の他、非破壊検査などに適用することも可能である。
 さらに、本例の信号処理方法により推定した患部260などの対象部分の位置,形状及び大きさなどに関する情報を入力信号として、患部260の位置を予測してもよい。
 1 制御部
 2 信号処理部
 3 放射線照射制御部
 4 透視画像取得部
 5 観測ベクトル生成部
 6 変数ベクトル更新部
 7 変数ベクトル計算部
 8 推定部
 90 放射線照射装置
 100 クライストロン
 110 冷却装置
 120 導波管
 130 真空ポンプ
 140 電子銃
 150 加速管
 160 偏向電磁石
 170 ターゲット
 180 平坦化フィルタ(散乱フォイル)
 190 モニタ線量計
 200 JAW(照射野形成部)
 210 MLC(マルチリーフコリメータ)
 220 OBI
 230 ベッド
 240 EPID
 250 パルスモジュレータ
 260 患部
 305 計測用放射線照射部 
  306 センサ部

Claims (16)

  1.  対象部分と前記対象部分とは異なる背景部分との透過的な重ね合わせにより表現される透視画像を処理する信号処理装置であって、
     複数の時刻について、前記対象部分を含む前記透視画像を取得する透視画像取得部と、
     前記複数の時刻のうちのある時刻tにおける前記対象部分に起因する前記透視画像の成分及び前記背景部分に起因する前記透視画像の成分を、前記複数の時刻における前記対象部分及び前記背景部分の分布情報の推定値と関連づけて計算し、前記複数の時刻における前記透視画像との整合性を評価する計算部と、
     前記計算部の評価結果に基づいて、前記時刻tにおける前記対象部分又は前記背景部分の分布情報の推定値を更新する更新部と、
    を備えることを特徴とする信号処理装置。
  2.  前記時刻tを他の時刻にずらしながら、前記計算部及び前記更新部における処理を行なうことにより、前記複数の時刻における前記対象部分又は前記背景部分の前記分布情報の推定値を更新する、
    ことを特徴とする、請求項1記載の信号処理装置。
  3.  前記分布情報には、前記対象部分及び前記背景部分の位置を示す位置情報及び輝度を示す輝度情報が含まれ、
     前記更新部は、前記位置情報又は前記輝度情報の推定値を更新する、
    ことを特徴とする、請求項1記載の信号処理装置。
  4.  前記複数の時刻における前記透視画像に基づいて、これらの時刻における各透視画像中の輝度情報を表す観測ベクトルbを生成する観測ベクトル生成部と、
     前記時刻tの透視画像における前記対象部分に起因する輝度情報を含む対象変数ベクトルIと、前記時刻tの透視画像における前記背景部分に起因する輝度情報を含む背景変数ベクトルIとを設定又は更新する変数ベクトル更新部と、
     前記複数の時刻における前記対象変数ベクトルI及び前記背景変数ベクトルIの前記透視画像中の位置情報を表現する関数Aを設定又は更新する関数更新部と、を備え、
     前記計算部は、前記観測ベクトルbと、前記対象変数ベクトルI及び前記背景変数ベクトルIを含む変数ベクトルIと、前記関数Aとの整合性を評価し、
     前記更新部は、前記計算部の評価結果に基づいて、前記関数A又は前記変数ベクトルIを更新する、
    ことを特徴とする、請求項3記載の信号処理装置。
  5.  前記計算部は、前記対象部分と前記背景部分との統計的な解析を行ない、
     前記更新部は、前記計算部における前記評価結果と前記統計的な解析結果とに基づいて、前記推定値を更新する、
    ことを特徴とする、請求項1記載の信号処理装置。
  6.  L(Lは自然数)個の対象部分と前記対象部分以外の背景部分との透過的な重ね合わせにより表現される透視画像を処理する信号処理装置であって、
     ある時刻tにおいて、前記対象部分に関し透視画像を取得する透視画像取得部と、
     上記透視画像取得部で取得された上記透視画像のM×N(M,Nは自然数)個の領域に着目し、前記時刻tと前記時刻tよりも過去のK(Kは自然数)個の時刻t1,t2,…,tKとにおいて取得した(K+1)個の透視画像の(K+1)×M×N個の前記着目する領域が有する輝度値からなる(K+1)時刻分の観測ベクトルb(t,t1,t2,…,tK)を生成する観測ベクトル生成部と、
     前記着目する領域中の前記L個の対象部分それぞれに起因する輝度値からなる変数ベクトル(以下、対象変数ベクトルIa1,Ia2,…,IaLという)と前記着目する領域中の前記背景部分が有する輝度値からなる変数ベクトル(以下、背景変数ベクトルIという)とを、前記観測ベクトルbと前記時刻tよりも過去のある時刻tpにおける透視画像の前記着目する領域に関するL個の変数ベクトルIa1(tp),Ia2(tp),…,IaL(tp)とからなる対象変数ベクトルIa(tp)及び背景変数ベクトルIb(tp)を用いてそれぞれ更新する変数ベクトル更新部と、
     上記観測ベクトルbと、前記変数ベクトル更新部により更新された前記対象変数ベクトルI及び背景変数ベクトルIからなる変数ベクトルIと、前記変数ベクトルIを用いて推定された前記時刻tにおける前記L個の対象部分及び背景部分の位置と前記K個の過去の時刻t1,t2,…,tKにおいて同様に推定されたそれぞれの時刻における前記L個の対象部分及び背景部分の位置とに関する(K+1)MN×(L+1)MNの係数行列Aとから、fを決定論的、もしくは確率論的な関数として、式
    Figure JPOXMLDOC01-appb-M000001

    により定義される評価値PIの評価を高めるように前記変数ベクトルIを計算する変数ベクトル計算部と、をそなえ、
     上記透視画像取得部での前記透視画像の取得、上記観測ベクトル生成部での前記観測ベクトルbの生成、上記変数ベクトル更新部での前記変数ベクトルIの更新、及び、上記変数ベクトル計算部での前記変数ベクトルIの計算を、時刻をずらしながら行なうように構成された、
    ことを特徴とする、信号処理装置。
  7.  前記変数ベクトル更新部が、
     前記対象変数ベクトルIにある初期値を設定し、
     前記背景変数ベクトルIに、前記観測ベクトルbから前記初期値が設定された対象変数ベクトルIを減算した結果を設定する、
    ことを特徴とする、請求項6記載の信号処理装置。
  8.  前記変数ベクトル更新部が、
     前記変数ベクトル計算部により計算された結果を用いて、前記変数ベクトルIを更新する、
    ことを特徴とする、請求項6または7に記載の信号処理装置。
  9.  前記変数ベクトル更新部が、
     前記変数ベクトル計算部により計算された背景変数ベクトルIにより、前記背景変数ベクトルIを更新し、
     前記観測ベクトルbから、当該更新された背景変数ベクトルIと過去のある時刻の(L-1)個の変数ベクトルIai,(i=1,2,…,L, ただし、i≠jでjは0<j<(L+1)を満たすある自然数)とを減算した結果を設定することにより、(L-1)個の前記変数ベクトルIajを更新し、
     前記観測ベクトルbから、当該更新された背景変数ベクトルIと前記更新された(L-1)個の変数ベクトルIaj,(j=1,2,…,L, ただし、j≠JでJは0<J<(L+1)を満たすある自然数)とを減算した結果を設定することにより、前記で更新されていない1個の前記変数ベクトルIaJを更新する、
    ことを特徴とする、請求項8記載の信号処理装置。
  10.  前記変数ベクトル計算部が、
     前記初期値を元に設定された変数ベクトルIと前記更新された変数ベクトルIとに基づいて、前記係数行列Aを推定する、
    ことを特徴とする、請求項7~9のいずれか1項に記載の信号処理装置。
  11.  前記変数ベクトル計算部が、
     前記更新される前の変数ベクトルIと前記更新された後の変数ベクトルIとに基づいて、前記係数行列Aを推定する、
    ことを特徴とする、請求項7~9のいずれか1項に記載の信号処理装置。
  12.  前記変数ベクトル計算部が、
     勾配法を用いて、前記変数ベクトルIを計算する、
    ことを特徴とする、請求項10または11に記載の信号処理装置。
  13.  前記変数ベクトル更新部によりある初期値が設定された対象変数ベクトルIと前記変数ベクトル更新部により更新された対象変数ベクトルIとに基づいて、前記対象部分の時間的な位置変動を推定する推定部をそなえる、
    ことを特徴とする、請求項6~12のいずれか1項に記載の信号処理装置。
  14.  前記変数ベクトル更新部により更新される前の対象変数ベクトルIと前記変数ベクトル更新部により更新された後の対象変数ベクトルIとに基づいて、前記対象部分の時間的な位置変動を推定する推定部をそなえる、
    ことを特徴とする、請求項6~12のいずれか1項に記載の信号処理装置。
  15.  対象部分と前記対象部分とは異なる背景部分との透過的な重ね合わせにより表現される透視画像を処理するコンピュータに対し、
     複数の時刻について、前記対象部分を含む前記透視画像を取得する透視画像取得機能と、
     前記複数の時刻のうちのある時刻tにおける前記対象部分に起因する前記透視画像の成分及び前記背景部分に起因する前記透過画像の成分を、前記複数の時刻における前記対象部分及び前記背景部分の分布情報の推定値と関連づけて計算し、前記複数の時刻における前記透視画像との整合性を評価する計算機能と、
     前記計算機能による評価結果に基づいて、前記時刻tにおける前記対象部分又は前記背景部分の分布情報の推定値を更新する更新機能と、
    を実現させるための信号処理プログラム。
  16.  対象部分と前記対象部分とは異なる背景部分との透過的な重ね合わせにより表現される透視画像を処理するコンピュータに、
     複数の時刻について、前記対象部分を含む前記透視画像を取得する透視画像取得機能と、
     前記複数の時刻のうちのある時刻tにおける前記対象部分に起因する前記透視画像の成分及び前記背景部分に起因する前記透過画像の成分を、前記複数の時刻における前記対象部分及び前記背景部分の分布情報の推定値と関連づけて計算し、前記複数の時刻における前記透視画像との整合性を評価する計算機能と、
     前記計算機能による評価結果に基づいて、前記時刻tにおける前記対象部分又は前記背景部分の分布情報の推定値を更新する更新機能と、
    を実現させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体。
PCT/JP2011/066007 2010-07-14 2011-07-13 信号処理装置、信号処理プログラム及び信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体 WO2012008500A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2012524581A JP5797197B2 (ja) 2010-07-14 2011-07-13 信号処理装置、信号処理プログラム及び信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体
EP11806830.3A EP2604187A4 (en) 2010-07-14 2011-07-13 SIGNAL PROCESSING DEVICE, SIGNAL PROCESSING PROGRAM, AND COMPUTER-READABLE RECORDING MEDIUM HAVING SIGNAL PROCESSING PROGRAM RECORDED THEREON
US13/740,859 US8837863B2 (en) 2010-07-14 2013-01-14 Signal-processing device and computer-readable recording medium with signal-processing program recorded thereon

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2010160009 2010-07-14
JP2010-160009 2010-07-14

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US13/740,859 Continuation US8837863B2 (en) 2010-07-14 2013-01-14 Signal-processing device and computer-readable recording medium with signal-processing program recorded thereon

Publications (1)

Publication Number Publication Date
WO2012008500A1 true WO2012008500A1 (ja) 2012-01-19

Family

ID=45469497

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2011/066007 WO2012008500A1 (ja) 2010-07-14 2011-07-13 信号処理装置、信号処理プログラム及び信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体

Country Status (4)

Country Link
US (1) US8837863B2 (ja)
EP (1) EP2604187A4 (ja)
JP (1) JP5797197B2 (ja)
WO (1) WO2012008500A1 (ja)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101537775B1 (ko) * 2014-06-25 2015-07-23 고신대학교 산학협력단 선형가속기의 정도 관리를 위한 템플릿 플래이트 및 정도 관리 템플릿 플래이트를 이용한 선형가속기의 정도 관리 방법
WO2016157457A1 (ja) * 2015-03-31 2016-10-06 国立大学法人東北大学 画像処理装置、画像処理方法、及び、画像処理プログラム
JP2017196009A (ja) * 2016-04-26 2017-11-02 コニカミノルタ株式会社 放射線撮影装置及び放射線撮影システム
JP6984908B2 (ja) * 2017-03-03 2021-12-22 国立大学法人 筑波大学 対象追跡装置
KR101876349B1 (ko) * 2017-04-20 2018-07-09 사회복지법인 삼성생명공익재단 치료기 제어 시스템 및 방법
US11369330B2 (en) * 2017-05-16 2022-06-28 Canon Medical Systems Corporation X-ray diagnosis apparatus
CN107918937B (zh) * 2017-12-06 2021-07-30 电子科技大学 一种基于光谱辐射的目标与背景的物理叠合方法
JP7147878B2 (ja) * 2019-01-21 2022-10-05 株式会社島津製作所 画像処理装置
JP7311319B2 (ja) * 2019-06-19 2023-07-19 ファナック株式会社 時系列データ表示装置
US12109059B2 (en) * 2020-05-19 2024-10-08 Konica Minolta, Inc. Dynamic analysis system, correction apparatus, storage medium, and dynamic imaging apparatus

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007105196A (ja) * 2005-10-13 2007-04-26 Fujifilm Corp 画像処理装置、画像処理方法およびそのプログラム
JP2008000190A (ja) * 2006-06-20 2008-01-10 Toshiba Corp X線診断装置およびx線診断装置におけるデータ処理方法
JP2008000456A (ja) * 2006-06-23 2008-01-10 Mitsubishi Heavy Ind Ltd 放射線治療装置制御装置および放射線照射方法
JP2008514352A (ja) * 2004-09-30 2008-05-08 アキュレイ インコーポレイテッド 運動中の標的の動的追跡
JP2010194053A (ja) * 2009-02-24 2010-09-09 Mitsubishi Heavy Ind Ltd 放射線治療装置制御装置および目的部位位置計測方法
JP2011018269A (ja) * 2009-07-10 2011-01-27 Nippon Telegr & Teleph Corp <Ntt> 半透明物体の動き検出装置および方法

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2632501B2 (ja) * 1994-07-12 1997-07-23 株式会社グッドマン 内臓器官の複数画像フレームを使用した生体内における内臓器官輪郭抽出方法
US6757423B1 (en) * 1999-02-19 2004-06-29 Barnes-Jewish Hospital Methods of processing tagged MRI data indicative of tissue motion including 4-D LV tissue tracking
JP3725418B2 (ja) * 2000-11-01 2005-12-14 インターナショナル・ビジネス・マシーンズ・コーポレーション 複数信号が混合される画像データから多次元信号を復元する信号分離方法、画像処理装置および記憶媒体
US7421101B2 (en) * 2003-10-02 2008-09-02 Siemens Medical Solutions Usa, Inc. System and method for local deformable motion analysis
US20050251029A1 (en) * 2004-04-21 2005-11-10 Ali Khamene Radiation therapy treatment plan
JP4064952B2 (ja) 2004-08-12 2008-03-19 三菱重工業株式会社 放射線治療装置および放射線治療装置の動作方法
US7555151B2 (en) * 2004-09-02 2009-06-30 Siemens Medical Solutions Usa, Inc. System and method for tracking anatomical structures in three dimensional images
CA2616306A1 (en) * 2005-07-22 2007-02-01 Tomotherapy Incorporated Method and system for processing data relating to a radiation therapy treatment plan
US7689021B2 (en) * 2005-08-30 2010-03-30 University Of Maryland, Baltimore Segmentation of regions in measurements of a body based on a deformable model
US7620205B2 (en) * 2005-08-31 2009-11-17 Siemens Medical Solutions Usa, Inc. Method for characterizing shape, appearance and motion of an object that is being tracked
US20070086639A1 (en) 2005-10-13 2007-04-19 Fujifilm Corporation Apparatus, method, and program for image processing
US7623679B2 (en) * 2006-12-13 2009-11-24 Accuray Incorporated Temporal smoothing of a deformation model
CN101909525B (zh) * 2008-01-11 2013-06-05 株式会社岛津制作所 图像处理方法、装置以及断层摄影装置
US8265363B2 (en) * 2009-02-04 2012-09-11 General Electric Company Method and apparatus for automatically identifying image views in a 3D dataset

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008514352A (ja) * 2004-09-30 2008-05-08 アキュレイ インコーポレイテッド 運動中の標的の動的追跡
JP2007105196A (ja) * 2005-10-13 2007-04-26 Fujifilm Corp 画像処理装置、画像処理方法およびそのプログラム
JP2008000190A (ja) * 2006-06-20 2008-01-10 Toshiba Corp X線診断装置およびx線診断装置におけるデータ処理方法
JP2008000456A (ja) * 2006-06-23 2008-01-10 Mitsubishi Heavy Ind Ltd 放射線治療装置制御装置および放射線照射方法
JP2010194053A (ja) * 2009-02-24 2010-09-09 Mitsubishi Heavy Ind Ltd 放射線治療装置制御装置および目的部位位置計測方法
JP2011018269A (ja) * 2009-07-10 2011-01-27 Nippon Telegr & Teleph Corp <Ntt> 半透明物体の動き検出装置および方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
See also references of EP2604187A4 *

Also Published As

Publication number Publication date
US8837863B2 (en) 2014-09-16
EP2604187A1 (en) 2013-06-19
JPWO2012008500A1 (ja) 2013-09-09
EP2604187A4 (en) 2016-09-07
JP5797197B2 (ja) 2015-10-21
US20130129255A1 (en) 2013-05-23

Similar Documents

Publication Publication Date Title
JP5797197B2 (ja) 信号処理装置、信号処理プログラム及び信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体
KR102070714B1 (ko) 오브젝트 위치 결정 장치, 오브젝트 위치 결정 방법, 오브젝트 위치 결정 프로그램, 및 방사선 치료 시스템
EP4095797B1 (en) Autonomous segmentation of three-dimensional nervous system structures from medical images
CN109074639B (zh) 医学成像系统中的图像配准系统和方法
JP5491174B2 (ja) 画像誘導型放射線治療のための画像の変形可能なレジストレーション
US9919164B2 (en) Apparatus, method, and program for processing medical image, and radiotherapy apparatus
CN111699021B (zh) 对身体中的目标的三维跟踪
EP3468668B1 (en) Soft tissue tracking using physiologic volume rendering
JP7113447B2 (ja) 医用画像処理装置、治療システム、および医用画像処理プログラム
JP5495886B2 (ja) 患者位置決めシステム
EP3629931B1 (en) Heatmap and atlas
EP4016470A1 (en) 3d morhphological or anatomical landmark detection method and device using deep reinforcement learning
US11645767B2 (en) Capturing a misalignment
JP7273416B2 (ja) 画像処理装置,画像処理方法および画像処理プログラム
Gill et al. Group-wise registration of ultrasound to CT images of human vertebrae
Siddique Towards Active Image Guidance in X-ray Fluoroscopy-guided Radiotherapeutic and Surgical Interventions

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 11806830

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2012524581

Country of ref document: JP

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 2011806830

Country of ref document: EP