[go: up one dir, main page]

WO2019230738A1 - コンピュータプログラム、画像処理装置、および画像処理方法 - Google Patents

コンピュータプログラム、画像処理装置、および画像処理方法 Download PDF

Info

Publication number
WO2019230738A1
WO2019230738A1 PCT/JP2019/021147 JP2019021147W WO2019230738A1 WO 2019230738 A1 WO2019230738 A1 WO 2019230738A1 JP 2019021147 W JP2019021147 W JP 2019021147W WO 2019230738 A1 WO2019230738 A1 WO 2019230738A1
Authority
WO
WIPO (PCT)
Prior art keywords
value
image
predetermined
image processing
images
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/JP2019/021147
Other languages
English (en)
French (fr)
Inventor
輝仁 城戸
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ehime University NUC
Original Assignee
Ehime University NUC
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 Ehime University NUC filed Critical Ehime University NUC
Priority to CN201980035304.2A priority Critical patent/CN112165901B/zh
Priority to EP19812179.0A priority patent/EP3811869B1/en
Priority to JP2020522222A priority patent/JP7236747B2/ja
Priority to US17/058,694 priority patent/US11813106B2/en
Publication of WO2019230738A1 publication Critical patent/WO2019230738A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
    • A61B6/5217—Devices 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
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A61B6/032—Transmission computed tomography [CT]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5205—Devices using data or image processing specially adapted for radiation diagnosis involving processing of raw data to produce diagnostic data
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/54—Control of apparatus or devices for radiation diagnosis
    • A61B6/541—Control of apparatus or devices for radiation diagnosis involving acquisition triggered by a physiological signal
    • G06T12/30—
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/48—Diagnostic techniques
    • A61B6/481—Diagnostic techniques involving the use of contrast agents
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/48—Diagnostic techniques
    • A61B6/486—Diagnostic techniques involving generating temporal series of image data
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
    • A61B6/503—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for diagnosis of the heart
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
    • A61B6/507—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for determination of haemodynamic parameters, e.g. perfusion CT
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/41—Medical
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/412—Dynamic

Definitions

  • the present invention relates to an image processing technique, and more particularly to a technique for analyzing a time-series dynamic image.
  • a method is known in which blood flow in a heart is analyzed using an image obtained by imaging the subject's heart after administering a contrast medium to the subject.
  • CTP Computed Tomography Perfusion
  • This CTP is a method for evaluating the state of myocardial blood flow, for example, by observing the first path of a contrast medium to the myocardium as a contrast effect.
  • TDC Time density curve
  • a CT image obtained by Dynamic imaging is analyzed, an input function and an output function are derived, and the derived input function and output are derived. Based on the function, an arrival time (hereinafter referred to as “AT”) that is the time at which the contrast agent arrives in a predetermined myocardial region and a base value that is a CT value that is a base of the pixel at the pixel position in the organ are calculated.
  • AT arrival time
  • an object of the present invention is to provide an analysis method that can analyze a time-series image by a simpler method than conventional methods.
  • Another object of the present invention is to provide an analysis method that reduces the exposure dose to the subject as compared with the prior art and ensures objectivity and quantitativeness in the analysis of time-series images.
  • the present invention has been made to solve at least a part of the problems described above, and can be realized as the following application examples.
  • reference numerals in parentheses and supplementary explanations in this column show correspondence with embodiments described later in order to help understanding of the present invention, and do not limit the present invention at all. .
  • the computer program according to Application Example 1 of the present invention is an image of a subject's organ after the contrast medium is administered, and is image data composed of a predetermined number of frames of CT (Computed Tomography) images in time series.
  • the first step (S350, S540) for specifying the change over time of the CT, and the inclination of the CT value with respect to the predetermined time and the predetermined inclination are specified based on the change over time of the CT value specified in the first step
  • the CT value based on the second step (S1100) and the predetermined slope specified in the second step A third step of approximating temporal changes in a predetermined function (S1300), and summarized in that for the execution.
  • the computer program of application example 2 in the present invention is the computer program of application example 1, wherein the predetermined inclination is an inclination of the CT value with respect to the predetermined time and is an inclination larger than a predetermined value. To do.
  • the computer program of application example 3 in the present invention is the computer program of application example 1 or application example 2, wherein the storage means includes first image data composed of CT images of a first predetermined number of frames, and the first program. Second image data composed of CT images of a second predetermined number of frames which is smaller than the predetermined number, and stored in the image processing device in the first step based on the first image data
  • a fourth step (S480) for specifying a maximum timing that is a timing at which the CT value becomes a maximum value is further executed, and the based on the maximum timing specified in the fourth step
  • a computer program according to Application Example 4 of the present invention is the computer program according to Application Example 3, in which the image processing apparatus performs one pulsation of the organ based on the maximum timing specified in the fourth step.
  • the gist is that the first step, the second step, and the third step are executed based on the second image data including the CT image of the corresponding second predetermined number of frames.
  • An image processing apparatus images an organ of a subject after administration of a contrast agent, and stores image data composed of a predetermined number of CT (Computed Tomography) images in time series.
  • a storage means 110
  • a first means 123, 133 for specifying a temporal change in CT value based on image data comprising CT images of the predetermined number of frames, and the first means specified by the first means
  • a second means 125, 137 for specifying a predetermined inclination that is an inclination of the CT value with respect to a predetermined time, and based on the predetermined inclination specified by the second means
  • a third means 125, 137 for approximating the temporal change of the CT value with a predetermined function.
  • the image processing method images an organ of a subject after administration of a contrast agent, and stores image data composed of a predetermined number of frames of CT (Computed Tomography) images in time series.
  • An image processing method executed by an image processing apparatus (100) provided with a storage means (110), wherein the image processing apparatus performs CT value aging based on image data comprising CT images of the predetermined number of frames.
  • a first step S350, S540
  • an inclination of the CT value with respect to a predetermined time and a predetermined inclination is specified.
  • the temporal change of the CT value is approximated by a predetermined function.
  • the third step (S1300) that will be summarized in that the execution.
  • FIG. 1 is an overall configuration diagram of an image processing apparatus 100 according to a first embodiment which is one embodiment of the present invention. It is explanatory drawing which illustrated an example of CT image. It is a flowchart which shows an input function derivation process. It is a flowchart which shows an input function specific process. It is a flowchart which shows an output function derivation process. It is a flowchart which shows an output function specific process. (A) is explanatory drawing which showed an example of the data structure of the input function table 700 of the input function data storage part 129, (B) shows an example of the data structure of the output function table 750 of the output function data storage part 139. FIG.
  • FIG. 1 is explanatory drawing which showed an example of 1st TDC of the aorta in Example 1
  • B is explanatory drawing which showed an example of 2nd TDC of the intramyocardial coronary artery in Example 1.
  • FIG. (A) is explanatory drawing which showed an example of 1st TDC of the aorta in Example 2
  • (B) is explanatory drawing which showed an example of 2nd TDC of the intramyocardial coronary artery in Example 2.
  • FIG. 1 is an overall configuration diagram of the image processing apparatus 100.
  • the image processing apparatus 100 is configured by, for example, a general-purpose computer system, and a computer program in which individual components or functions in the image processing apparatus 100 described below are stored in a computer-readable recording medium or the like. It is realized by executing.
  • the image processing apparatus 100 according to the present embodiment includes an image data storage unit 110, an input function derivation unit 120, an input function data storage unit 129, an output function derivation unit 130, an output function data storage unit 139, and a blood flow analysis. And a processing unit 150.
  • an input device 160 and a display device 170 are connected to the image processing apparatus 100.
  • the image processing apparatus 100 is an apparatus that performs quantitative analysis of blood flow of an organ using a CT image of the organ of the subject imaged after administering a contrast agent to the subject.
  • the image processing apparatus 100 targets a plurality of CT frame images of the same phase obtained by imaging the subject's heart in synchronization with the electrocardiogram, and analyzes the heart CT images in units of pixels.
  • An input function, an output function, and the like used for calculating a quantitative value of the blood flow in the heart and in the vicinity of the heart are specified.
  • the image processing apparatus 100 estimates the blood flow rate of the heart or the like based on the change in the pixel value of the CT image by the contrast agent administered to the vein, and the contrast agent reaches the target site, for example.
  • the pixel value before the pixel value starts to change is set as the base value, or the time when the pixel value starts to change is set as AT.
  • this contrast agent iopromide shall be utilized, for example.
  • the image data storage unit 110 is a general memory in a computer system.
  • the image data storage unit 110 stores image data of a heart and a CT image near the heart obtained by imaging the heart.
  • a CT image is a CT image of a plurality of time-series frames in which an organ of a subject after administration of a contrast agent is imaged, and is a three-dimensional image (three-dimensional image) composed of a plurality of slice images (short-axis tomographic images). Voxel data), details of which will be described later.
  • the input function deriving unit 120 Based on the image data stored in the image data storage unit 110, the input function deriving unit 120 affects the contrast agent flowing into the organ of the subject on the pixel value (hereinafter referred to as “CT value”) of the CT image.
  • CT value pixel value
  • the input function related to the change is derived.
  • the input function deriving unit 120 includes an ROI setting unit 121, an ROI temporal change specifying unit 123, and an input function conversion processing unit 125 described below.
  • the input function data storage unit 129 stores data related to the input function derived by the input function deriving unit 120 as an input function table 700 described later.
  • the ROI setting unit 121 sets a region of interest (ROI), which is an area for specifying an input function in a CT image, in accordance with an operation via an input device 160 by an analyst such as a doctor. Specifically, the ROI setting unit 121 reads the image data from the image data storage unit 110 and causes the display device 170 to display a frame image of a slice in which the target region is reflected, and the analyst uses the input device 3. The ROI is set on the image displayed on the display device 170 in accordance with the input made in the above. Note that this ROI position is common to all frames of the same slice.
  • ROI region of interest
  • the ROI temporal change specifying unit 123 specifies temporal changes in CT values in the ROI based on the image data of CT images of a plurality of time-series frames. Specifically, the ROI temporal change specifying unit 123 generates a TDC of the ROI value determined from the CT value in the ROI as the temporal change of the CT value of a predetermined region in the CT image.
  • the ROI value is a value (statistical value) obtained by processing the CT value in the ROI with a statistical algorithm. For example, any one of an average value, a mode value, a median value, a maximum value, and a minimum value is taken. To get.
  • the input function processing unit 125 Based on the TDC generated by the ROI temporal change specifying unit 123, the input function processing unit 125 reaches the region where the MT and ROI are set at which the CT value in the ROI becomes the maximum value Max.
  • the AT that is the time when the image is taken, the MB that is the imaging section of the imaging period for obtaining the second TDC described later, the input function, and the like are specified.
  • the input function processing unit 125 specifies an input function and the like after performing smoothing processing so that the TDC becomes a smooth curve in the TDC generated by the ROI temporal change specifying unit 123. Further, the input function conversion processing unit 125 also specifies an upper limit frame Fa that is a predetermined frame after the CT value has rapidly increased, and BL that is a base value described later.
  • the output function deriving unit 130 is an output function related to a change in the CT value of the CT image due to the contrast agent flowing into the blood vessel of the organ of the subject. Is derived. Note that the output function deriving unit 130 derives an output function in units of pixels instead of in units of ROIs, and includes a target pixel extracting unit 131, a CT temporal change specifying unit 133, and an output functioning processing unit 137. doing.
  • the output function data storage unit 139 stores data relating to the output function derived by the output function deriving unit 130 as an output function table 750 described later.
  • the target pixel extraction unit 131 extracts a pixel from which an output function is derived, and specifies its position on the frame image. Specifically, the target pixel extraction unit 131 selects pixels whose CT values are within a predetermined range (for example, 50 to 160) in units of slices in all time-series CT images of a plurality of frames as extraction processing. . In addition, in this extraction process, when a region that is not a small target is included in the region of the target pixel of interest, according to the size of the region that is not the target (for example, the number of pixels is a predetermined number or less) Therefore, it is possible to execute processing for converting this region into a target pixel by regarding this region as a processing defect.
  • a predetermined range for example, 50 to 160
  • the target pixel extraction unit 131 identifies the position of the pixel extracted for each slice as the target pixel position (intra-organ pixel position) that is the target for calculating the output function.
  • the CT temporal change specifying unit 133 specifies the temporal change of the CT value of the pixel at the target pixel position (pixel position in the organ) based on the image data of the CT images of a plurality of time-series frames. Specifically, the CT temporal change specifying unit 133 generates a TDC as the temporal change of the CT value of the target pixel position for each slice.
  • the output function conversion processing unit 137 specifies an output function or the like based on the TDC generated by the CT temporal change specifying unit 133. Specifically, the output functioning processing unit 137 calculates the slope value between the CT values based on the TDC that is the temporal change specified by the CT temporal change specifying unit 133, and among the calculated slope values, The output function is specified on the basis of the slope of a predetermined value or more.
  • the blood flow analysis processing unit 150 performs a quantitative analysis of the blood flow based on the input function and the output function. Specifically, the blood flow analysis processing unit 150 performs a predetermined analysis process with reference to the input function table 700 of the input function data storage unit 129 and the output function table 750 of the output function data storage unit 139.
  • a technique for quantitative analysis of the blood flow for example, using the Patlake Plot method or the Deconvolution method, MBF (Myocardial Blood Flow: Myocardial Blood Volume), MBV (Myocardial Blood Volume), MTT (Mean Blood) (Transit Time: average transit time) and the like.
  • a predetermined analysis process is executed based on the slope a or function F data of the function F in the input function table 700, which will be described later, and the slope ⁇ or data of the function F in the output function table 750.
  • the result of quantitative analysis by the blood flow analysis processing unit 150 is displayed on the display device 170.
  • the result of quantitative analysis is displayed side by side or overlaid with the CT image, and the two images are compared, thereby evaluating coronary artery stenosis and myocardial ischemia in the case of myocardial infarction or angina. It can be done at the same time.
  • the identified MBF when the identified MBF is compared with the normal MBF, and the identified MBF is lower than the normal MBF, it can be determined that the ischemic state is present, and consequently, in the coronary artery of the myocardium. Evaluation of coronary artery stenosis and myocardial ischemia during myocardial infarction or angina can be performed.
  • FIG. 2 is an explanatory diagram illustrating an example of a CT image.
  • the CT image shown in FIG. 2 is a frame image 210 in which the heart is imaged in synchronization with the electrocardiogram and a predetermined number (for example, 30) of heartbeats are in phase.
  • One frame is composed of a plurality of slice images, and one frame of data is composed of three-dimensional image data.
  • FIG. 2 illustrates a frame image of one slice.
  • the CT image may be one in which imaging is started immediately after the contrast medium is administered to the subject, and may include the first frame to the 30th frame in the order of imaging.
  • FIG. 3 is a flowchart showing the input function derivation process.
  • FIG. 4 is a flowchart showing the input function specifying process.
  • the ROI setting unit 120 is an image of a frame image of a slice image in which a specific area is captured from the image data storage unit 110 in accordance with the selection process by the analyst.
  • Data is read (S310), and when a frame image in which a specific area is clearly shown in the selected slice image is displayed on the display device 5, an ROI is set according to the operation of the analyst (S330).
  • the ROI temporal change specifying unit 123 calculates an ROI value for each frame image for which the ROI is set, and draws a TDC based on the ROI value (S350).
  • the input function conversion processing unit 125 executes an input function specifying process described later (S390).
  • the input function processing unit 125 defines a range of frames from which an input function is derived based on the TDC generated by the ROI temporal change specifying unit 123.
  • the upper limit frame Fa to be identified is specified (S405), the variable n indicating the frame number is set to 1 (S410), and 1 is added to n (S420).
  • the input function conversion processing unit 125 approximates the ROI value of the frame before the nth frame to one straight line L (S430), and secondarily calculates the ROI value of the frame after the nth frame and before the upper limit frame Fa.
  • the function F is approximated (S440).
  • the input function conversion processing unit 125 calculates the least square error of the ROI value of the frame ⁇ before the straight line L and the nth frame, and the second function F and the frames after the nth frame and before the upper limit frame Fa.
  • the least square error with the ROI value is calculated, and the sum of these errors is calculated (S450).
  • the ROI value (height of the straight line L (Y intercept)) in the frame designated as AT is BL
  • the period corresponding to a predetermined heart rate before and after MT is the imaging section MB
  • the input function and the slope of the input function a is specified.
  • FIG. 5 is a flowchart showing the output function derivation process.
  • FIG. 6 is a flowchart showing the output function specifying process.
  • the target pixel extraction unit 131 selects one slice image, reads the image data of the frame image of the slice from the image data storage unit 110 (S510), Among all the pixels of the selected slice, pixels whose CT values in all frame images satisfy the above-mentioned conditions are regarded as pixels in the specific area, and these pixels are extracted as target pixels (S520), and the extracted target One pixel is selected from the pixels (S530).
  • the CT temporal change specifying unit 133 draws a TDC for the pixel selected by the target pixel extracting unit 131 (S540).
  • the output function conversion processing unit 137 specifies an output function or the like by executing an output function specifying process described later based on the CT value (S550). Then, the output function deriving unit 130 performs the processing from steps S530 to S550 for all target pixels (S580) and performs the processing from steps S510 to S580 for all slices (S590).
  • the output function conversion processing unit 137 calculates the slope ⁇ between the plots in the TDC generated by the CT temporal change specifying unit 133 (S610). Only those slopes that are greater than or equal to a predetermined value (for example, greater than or equal to 0.2) are extracted (S620). Then, the output function conversion processing unit 137 approximates with a predetermined function based on the extracted inclination ⁇ (S630), and specifies an output function (S650).
  • a predetermined value for example, greater than or equal to 0.2
  • FIGS. 7A and 7B are explanatory diagrams showing examples of data structures of the input function table 700 and the output function table 750.
  • FIG. 7A shows data in the input function table 700 of the input function data storage unit 129. It is explanatory drawing which showed an example of a structure, (B) is explanatory drawing which showed an example of the data structure of the output function table 750 of the output function data storage part 139.
  • the input function table 700 includes AT711, MT712, BL713, function function slope a714, functions L and F715, ROI value 717, and MB718 as data items. It is out.
  • the blood flow analysis processing unit 150 inputs an input obtained by correcting the ROI value 717 or the functions L and F715 so that the AT 711 and BL 715 become the origin. It may be used as a function.
  • the output function table 750 includes a slice number 761, a pixel position 762, a slope ⁇ 763, and a function F765 as data items.
  • the pixel position 762 is data indicating the positions of the two CT values that are the basis for calculating the inclination ⁇ 763
  • the function F765 is an expression of the function F that is an approximate expression calculated from the inclination ⁇ 763.
  • the contrast medium diluted to a predetermined concentration with physiological saline for example, the contrast medium diluted four times with physiological saline
  • imaging using an undiluted contrast agent is performed.
  • imaging is performed by the test injection method, whereby a first CT image that is a CT image of the heart including the aorta of the subject is captured.
  • the image processing apparatus 100 by executing the input function derivation process as described above by the image processing apparatus 100 based on the first CT image, the first TDC of the aorta is generated, and the MB, the input function, and the like are specified.
  • imaging using a non-diluted contrast agent is executed, whereby a second CT image, which is a CT image of the heart including the subject's aorta and intramyocardial coronary artery, is captured.
  • the output function deriving process as described above is executed by the image processing apparatus 100 based on the second CT image, thereby generating the second TDC of the intramyocardial coronary artery and specifying the output function and the like.
  • the image processing apparatus 100 performs a quantitative analysis of the blood flow rate in the intramyocardial coronary artery.
  • FIGS. 8A and 8B are explanatory diagrams illustrating an example of the TDC in the present embodiment
  • FIG. 8A is an explanatory diagram illustrating an example of the first TDC of the aorta in the present embodiment.
  • 8A shows the case where the ROI is set for the aorta in the first CT image in the vicinity of the aorta of the heart.
  • FIG. 8B shows the target pixel with a predetermined CT value (in this embodiment, 50 to 150 HU)).
  • the first TDC of the aorta is generated by the image processing of the image processing apparatus 100 based on the first CT image of the heart obtained by imaging by the test injection method.
  • MB etc. are specified based on 1 TDC.
  • the ROI temporal change specifying unit 123 of the input function deriving unit 120 first applies the ROI set by the ROI setting unit 121 to other frame images 210 of the same slice, and for each frame image, A ROI value representing the ROI is determined based on the CT value of the pixel in the ROI.
  • the ROI temporal change specifying unit 123 generates a first TDC by plotting ROI values determined for all frames of the same slice.
  • the ROI temporal change specifying unit 123 detects the maximum value Max of the CT value in the first TDC, or the first peak value after the curve has risen significantly (the CT value suddenly increases).
  • the ROI temporal change specifying unit 123 can set a fixed ratio (for example, 70%, 80%, 90%, etc.) from the minimum value Min of the ROI value in the TDC to the maximum value Max or the peak value. In this embodiment, the value is 70%).
  • the input function processing unit 125 is a frame before the frame in which the maximum value Max or the peak value is detected, and takes a certain ratio value specified by the ROI temporal change specifying unit 123 as described above. Let the frame be the upper limit frame Fa. Next, the input function conversion processing unit 125 performs linear approximation by applying the least square method or the like to the ROI value before the nth (n is from 2 to the upper limit frame Fa-1) frame Fn, and the straight line L In addition, the least square method or the like is applied to the ROI value from the nth frame Fn to the upper limit frame Fa to approximate a quadratic function to derive the quadratic function F.
  • the input function processing unit 125 calculates the sum of squares of errors between the straight line L and the ROI values before the nth frame Fn, and calculates the error between the quadratic function F and the ROI values after the nth frame Fn. Calculate the sum of squares (residual sum of squares). Then, the input function conversion processing unit 25 performs the above-described processing for all n, specifies n that minimizes the sum of the squares of these errors, and straight lines L and 2 in the case of n that minimizes this error. The next function F is specified as an approximated input function, and the slope a of the input function is specified.
  • the frame in which the ROI value is the maximum value Max is MT
  • the frame that is the intersection (boundary) of the straight line L and the quadratic function F is AT
  • the ROI value is BL
  • the section corresponding to 3 heartbeats before and after MT is MB.
  • the input function conversion processing unit 125 identifies MT, AT, BL, MB, slope a, and input function of the aortic region based on the first TDC.
  • the height of the straight line L (Y intercept) may be set to BL.
  • MB is an imaging section that is a timing for performing imaging for obtaining a second CT image.
  • the imaging interval MB in this embodiment includes one heartbeat and MT before MT. The interval corresponds to the latter two heartbeats.
  • the second TDC of the intramyocardial coronary artery is generated by the image processing of the image processing apparatus 100 based on the second CT image obtained by the imaging in the MB, and the second TDC is included in the second TDC. Based on this, an output function or the like is specified.
  • blood mainly flows into the intramyocardial coronary artery in the diastole of the heart accompanying the heartbeat. Therefore, the CT value changes between the diastole and the non-diastolic phase of the heart, and the second TDC has a step difference. It becomes a plot that is formed.
  • the output function when specifying the output function based on the second TDC, it is possible to specify a highly accurate output function by approximation based on the value of the slope that is equal to or greater than a predetermined value. Specifically, the values of the slopes between the plots in the second TDC based on the imaging in the imaging section MB are calculated, and the values of the slopes ⁇ 1 to ⁇ 3 that are equal to or greater than a predetermined value among all the slope values are extracted. The output function is specified by the output function specifying process based on the inclinations ⁇ 1 to ⁇ 3.
  • the change in the slope value during one heartbeat in the second TDC along the relationship between the diastole and the non-diastolic phase of the heart has a small slope in the first half and a large slope in the second half. It is an aspect which becomes the value of.
  • the slope values used in the output function specifying process are slopes ⁇ 1 to ⁇ 3, which are large slope values in the latter half of one heartbeat. That is, in other words, the output function specifying process is a process using the slope value of the latter half of the slope during one heartbeat in the second TDC.
  • the output function may be specified using an average value of the gradients ⁇ 1 to ⁇ 3, or after specifying each approximate expression of the gradients ⁇ 1 to ⁇ 3, various identification methods may be used. An output function may be specified.
  • FIG. 9A and 9B are explanatory diagrams illustrating an example of TDC in the second embodiment
  • FIG. 9A is an explanatory diagram illustrating an example of the first TDC of the aorta in the second embodiment.
  • FIG. The configuration of the image processing apparatus 100 in the second embodiment is the same as that in the first embodiment, and the flow of the input function derivation process and the output function derivation process in the second embodiment is substantially the same as that in the first embodiment.
  • ⁇ Configuration of image processing apparatus 100>, ⁇ Contents of input function derivation process>, ⁇ Contents of output function derivation process>, and the like are omitted.
  • the contrast agent that is not diluted after the heart is imaged by the test injection method.
  • the MB specified based on the first TDC in the second embodiment is a section corresponding to one heartbeat, and the second TDC is also generated corresponding to one heartbeat.
  • the inclination ⁇ 4 in the second TDC based on the imaging in the imaging section MB is extracted, and the output function is specified by the output function specifying process based on the inclination ⁇ 4.
  • MB which is an imaging section in the second embodiment, is a section corresponding to one heartbeat after MT in order to obtain the CT value of the intramyocardial coronary artery in the frame when the contrast medium surely flows into the intramyocardial coronary artery.
  • the image data storage unit 110 includes the first image data including the first CT images of the first predetermined number of frames and the first predetermined number.
  • Second image data consisting of a second CT image of a second predetermined number of frames, which is smaller than that of the second CT image, and based on the MT specified from the first CT image obtained by the test injection method using a diluted contrast agent.
  • an MB that is an imaging section corresponding to a predetermined heart rate is specified, imaging using an undiluted contrast agent in MB is performed, and the image processing apparatus 100 is based on the obtained second CT image.
  • the image processing can be executed.
  • the image processing apparatus 100 uses the second CT image and the second TDC having the number of frames corresponding to one pulsation based on the MB that is the imaging section. Image processing can be performed. Therefore, it is possible to perform image processing with a CT image by the CTP examination in the imaging section MB, which is a period corresponding to a small heart rate, as compared with the conventional case of performing a CTP examination corresponding to a period of 20 to 30 heartbeats. Therefore, it is possible to reduce the exposure dose in the subject as compared with the conventional case and execute image processing based on a simple examination.
  • the section corresponding to three or one heartbeat before and after the MT based on the first CT image is the imaging section MB.
  • the present invention is not limited to this, and the MB is as follows. Also good. For example, a period corresponding to 2 heartbeats or 4 to 10 heartbeats before and after MT may be set as the imaging section MB. Even in such shooting based on MB or imaging timing, it is possible to reduce the exposure dose in the subject as compared with the conventional case and execute image processing based on a simple examination.
  • the present invention is not limited to this, and 0.2.
  • a numerical value other than may be used as the threshold value.
  • a numerical value based on a statistical value calculated from a result of a plurality of shootings may be used as the threshold value. With such a threshold value, it is possible to perform image processing with higher accuracy, and thus quantitative analysis of blood flow with higher accuracy.
  • the output function specifying process shown in FIG. 6 extracts only the slope that is 0.2 or more in the calculated slope, that is, the slope in the diastole of the heart (of the slope during one heartbeat).
  • the present invention is not limited to this, and the output function may be specified based on a change in CT value during the non-diastolic period of the heart. Specifically, for example, a slope that is less than 0.2 in the second TDC is extracted, and a difference based on a difference between changes in each CT value (change in CT value in the non-diastolic period of the heart) based on the extracted slope.
  • the output function may be specified by a method or the like.
  • the output function specifying process shown in FIG. 6 is used in the output function derivation process.
  • the present invention is not limited to this, and the output function shown in FIG. A process such as a specific process may be used.
  • the input function may be specified from the slope as in the output function specifying process.
  • the Dynamic imaging is performed by the test injection method using the contrast medium in which the contrast medium is diluted with physiological saline.
  • imaging may be performed by a bolus tracking method in which an ROI is set for a target region and a CT value (ROI value) of the ROI is monitored in real time.
  • Dynamic imaging may be performed by a test injection method using an undiluted contrast agent, and Dynamic imaging may be performed by a so-called test bolus tracking method that combines a test injection method and a bolus tracking method. May be. That is, in the above-described embodiment, in order to optimize the timing for performing imaging for capturing the second CT image, the first CT image can be captured using various techniques.
  • the first CT image is obtained by imaging with the test injection method using the diluted contrast agent, and then the second CT image is obtained by imaging using the undiluted contrast agent.
  • the present invention is not limited to this.
  • a technique may be employed in which a second CT image is captured by capturing a plurality of cardiac phases only at a predetermined heart rate during execution of imaging by the test injection method for capturing the first CT image.
  • the target pixel extraction unit 131 executes the extraction process as described above.
  • the present invention is not limited to this, and the following process may be executed for each pixel.
  • the target pixel extraction unit 131 may select a pixel based on the amount of change in the CT value of each pixel, for example, the difference between the maximum value and the minimum value, in a time-series CT image of a plurality of frames. .
  • the target pixel extraction unit 131 obtains a difference between the maximum value and the minimum value of the CT values in all frame images in one slice, and the difference is a predetermined value (for example, 50 to 150). Pixels may be extracted.
  • the target pixel extraction unit 131 may execute such extraction processing in addition to the above-described extraction processing, or may execute only one of them.
  • the upper limit frame Fa is specified by the above-described method in the input function specifying process, but the present invention is not limited to this.
  • the change rate of the first TDC may be obtained, and any frame between the frame having the maximum change rate and the frame having the change rate of 0 may be set as the upper limit frame Fa.
  • the first TDC is functionalized with two functions of a straight line and a quadratic curve, but the present invention is not limited to this.
  • the line before the nth frame Fn may be approximated by a straight line L
  • the distribution of ROI values or CT values after the nth frame Fn may be approximated by a straight line or a multi-order function of third or higher order.
  • the distribution of all ROI values or CT values up to the upper limit frame Fa may be approximated to a function represented by a multidimensional polynomial.
  • the first TDC may be approximated by three or more functions.
  • the second frame to the upper limit frame Fa-1 are divided into three or more sections, and each section is approximated by a predetermined function. Also good.
  • the ROI value or CT value may temporarily decrease immediately before the curve rises significantly (ROI value or CT value increases rapidly).
  • a section (first section) in which the ROI value or CT value is almost constant and linearly approximated a section in which the ROI value or CT value decreases (second section), and the subsequent ROI You may divide into the section (3rd section) where a value or CT value rises rapidly.
  • the first section may be approximated by a straight line
  • the second section may be approximated by a second-order or higher function
  • the third section may be approximated by another second-order or higher-order function.
  • the frame serving as the boundary between the first section and the second section or the frame serving as the boundary between the second section and the third section is defined as AT
  • the ROI value or CT of the frame serving as this AT The value may be a base value.
  • the first TDC may be approximated by one function.
  • the normal TMD or the cumulative distribution function may be fitted to the first TDC, and the first TDC may be approximated by a single function.
  • the standard deviation (SD) and average value of the normal distribution may be selected so as to best fit the rising curve of the first TDC.
  • the frame closest to ⁇ 3SD of the normal cumulative distribution function approximated to the first TDC may be set as AT, and the ROI value or CT value of the frame serving as the AT may be set as the base value.
  • the function specifying process approximates the smoothed first TDC with a function.
  • the present invention is not limited to this, and the first TDC can be obtained by a method other than fitting to an m-order function as described above. May be smoothed.
  • an average value with surrounding pixels may be obtained and smoothed using the average value, or the first TDC may be smoothed by a moving average.
  • a function approximation may be performed on the first TDC that has not been smoothed.
  • the result of quantitative analysis by the blood flow analysis processing unit 150 is displayed on the display device 170 as described above, but the present invention is not limited to this.
  • each segment may be displayed in a display mode corresponding to the blood flow rate.
  • the 3D image of the organ based on the image data preserve
  • the CT image may be displayed as a 3D image based on the coordinate information of each pixel.
  • DESCRIPTION OF SYMBOLS 100 ... Image processing apparatus, 110 ... Image data memory

Landscapes

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

Abstract

本発明は、従来に比べてより簡便な手法で時系列画像を解析できる解析手法を提供することを目的としている。本発明の適用例の一つであるコンピュータプログラムは、造影剤が投与された後の被験者の臓器を撮像したものであって、時系列に沿った所定数のフレームのＣＴ画像からなる画像データを記憶する記憶手段を備えた画像処理装置のためのコンピュータプログラムであって、画像処理装置に、所定数のフレームのＣＴ画像からなる画像データに基づいて、ＣＴ値の経時変化を特定する第１ステップと、第１ステップで特定されたＣＴ値の経時変化に基づいて、所定時間に対するＣＴ値の傾きであって所定の傾きを特定する第２ステップと、第２ステップで特定された所定の傾きに基づいて、ＣＴ値の経時変化を所定の関数で近似する第３ステップと、を実行させる。

Description

コンピュータプログラム、画像処理装置、および画像処理方法
　本発明は、画像処理技術に関し、特に時系列のダイナミック画像を解析する技術に関する。
　被験者に造影剤を投与した後、その被験者の心臓を撮像した画像を使って、心臓の血流を解析する手法が知られている。このような心筋血流の解析手法の一つとして、冠動脈ＣＴ検査時に追加して行われるＣＴＰ（Ｃｏｍｐｕｔｅｄ　Ｔｏｍｏｇｒａｐｈｙ　Ｐｅｒｆｕｓｉｏｎ）がある。このＣＴＰは、例えば、造影剤の心筋へのファーストパスを造影効果として観察することで、心筋血流の状態を評価する方法である。また、このＣＴＰにおける評価法としては、徐々に染まっていく心筋を複数心拍で撮像するＤｙｎａｍｉｃ撮影によって得られたＴＤＣ（Ｔｉｍｅ　ｄｅｎｓｉｔｙ　ｃｕｒｖｅ）を解析することで、心筋血流を定量評価するものがある。
　このようにＴＤＣの解析によって心筋血流等を定量評価する手法としては、例えば、Ｄｙｎａｍｉｃ撮影によって得られたＣＴ画像を解析して、入力関数および出力関数を導出して、導出した入力関数および出力関数に基づいて所定の心筋領域での造影剤が到達した時刻であるアライバルタイム（以下、「ＡＴ」とする。）や臓器内画素位置の画素のベースとなるＣＴ値であるベース値を算出したうえで、心筋血流を定量評価する方法がある（例えば、特許文献１を参照。）。
国際公開２０１６／００９９５７号
　しかしながら、上述した特許文献１に記載された方法では、少なくとも３０～４０心拍数におけるＤｙｎａｍｉｃ撮影によって得られた複数フレームのＣＴ画像を解析する必要があり、解析処理が煩雑であるため、より簡便な手法で時系列画像を解析できるような解析方法が望まれている。また、少なくとも３０～４０心拍数に応じた時間のＣＴＰ検査が必要であるため、被験者に対しての被ばく量が多くなってしまう、という問題があった。そのため、撮像タイミングを的確に把握し、撮像時間を短くして、被験者に対する被ばく量を少なくする、ということも望まれている。
　そこで、本発明の目的は、従来に比べてより簡便な手法で時系列画像を解析できる解析手法を提供することである。また、本発明の別の目的は、従来に比べて被験者に対する被ばく量を少なくするとともに、時系列画像の解析において客観性及び定量性を担保した解析手法を提供することである
　本発明は、上述の課題の少なくとも一部を解決するためになされたものであり、以下の適用例として実現することが可能である。なお、本欄における括弧内の参照符号や補足説明等は、本発明の理解を助けるために、後述する実施例との対応関係を示したものであって、本発明を何ら限定するものではない。
　本発明における適用例１のコンピュータプログラムは、造影剤が投与された後の被験者の臓器を撮像したものであって、時系列に沿った所定数のフレームのＣＴ（Ｃｏｍｐｕｔｅｄ　Ｔｏｍｏｇｒａｐｈｙ）画像からなる画像データを記憶する記憶手段（１１０）を備えた画像処理装置（１００）のためのコンピュータプログラムであって、前記画像処理装置に、前記所定数のフレームのＣＴ画像からなる画像データに基づいて、ＣＴ値の経時変化を特定する第１ステップ（Ｓ３５０、Ｓ５４０）と、前記第１ステップで特定された前記ＣＴ値の経時変化に基づいて、所定時間に対する前記ＣＴ値の傾きであって所定の傾きを特定する第２ステップ（Ｓ１１００）と、前記第２ステップで特定された前記所定の傾きに基づいて、前記ＣＴ値の経時変化を所定の関数で近似する第３ステップ（Ｓ１３００）と、を実行させることを要旨とする。
　本発明における適用例２のコンピュータプログラムは、適用例１のコンピュータプログラムであって、前記所定の傾きは、前記所定時間に対する前記ＣＴ値の傾きであって所定値より大きい傾きであることを要旨とする。
　本発明における適用例３のコンピュータプログラムは、適用例１または適用例２のコンピュータプログラムであって、　前記記憶手段は、第１所定数のフレームのＣＴ画像からなる第１画像データと、前記第１所定数に比べて少ない第２所定数のフレームのＣＴ画像からなる第２画像データと、を記憶しており、前記画像処理装置に、前記第１画像データに基づき前記第１ステップで特定された前記ＣＴ値の経時変化において、前記ＣＴ値が最大値となるタイミングである最大タイミングを特定する第４ステップ（Ｓ４８０）、をさらに実行させ、前記第４ステップで特定された前記最大タイミングに基づく前記第２所定数のフレームのＣＴ画像からなる前記第２画像データに基づいて、前記第１ステップ、前記第２ステップ、および前記第３ステップを実行させる、ことを要旨とする。
　本発明における適用例４のコンピュータプログラムは、適用例３のコンピュータプログラムであって、前記画像処理装置に、前記第４ステップで特定された前記最大タイミングに基づいた前記臓器の１回の拍動に相当する前記第２所定数のフレームのＣＴ画像からなる前記第２画像データに基づいて、前記第１ステップ、前記第２ステップ、および前記第３ステップを実行させる、ことを要旨とする。
　本発明における画像処理装置は、造影剤が投与された後の被験者の臓器を撮像したものであって、時系列に沿った所定数のフレームのＣＴ（Ｃｏｍｐｕｔｅｄ　Ｔｏｍｏｇｒａｐｈｙ）画像からなる画像データを記憶する記憶手段（１１０）と、前記所定数のフレームのＣＴ画像からなる画像データに基づいて、ＣＴ値の経時変化を特定する第１手段（１２３、１３３）と、前記第１手段によって特定された前記ＣＴ値の経時変化に基づいて、所定時間に対する前記ＣＴ値の傾きであって所定の傾きを特定する第２手段（１２５、１３７）と、前記第２手段によって特定された前記所定の傾きに基づいて、前記ＣＴ値の経時変化を所定の関数で近似する第３手段（１２５、１３７）と、を備えることを要旨とする。
　本発明における画像処理方法は、造影剤が投与された後の被験者の臓器を撮像したものであって、時系列に沿った所定数のフレームのＣＴ（Ｃｏｍｐｕｔｅｄ　Ｔｏｍｏｇｒａｐｈｙ）画像からなる画像データを記憶する記憶手段（１１０）を備えた画像処理装置（１００）が実行する画像処理方法であって、前記画像処理装置が、前記所定数のフレームのＣＴ画像からなる画像データに基づいて、ＣＴ値の経時変化を特定する第１ステップ（Ｓ３５０、Ｓ５４０）と、前記第１ステップで特定された前記ＣＴ値の経時変化に基づいて、所定時間に対する前記ＣＴ値の傾きであって所定の傾きを特定する第２ステップ（Ｓ１１００）と、前記第２ステップで特定された前記所定の傾きに基づいて、前記ＣＴ値の経時変化を所定の関数で近似する第３ステップ（Ｓ１３００）と、を実行することを要旨とする。
本発明の一つの実施例である実施例１における画像処理装置１００の全体構成図である。 ＣＴ画像の一例を例示した説明図である。 入力関数導出処理を示すフローチャートである。 入力関数特定処理を示すフローチャートである。 出力関数導出処理を示すフローチャートである。 出力関数特定処理を示すフローチャートである。 （Ａ）は入力関数データ記憶部１２９の入力関数テーブル７００のデータ構成の一例を示した説明図であり、（Ｂ）は出力関数データ記憶部１３９の出力関数テーブル７５０のデータ構成の一例を示した説明図である。 （Ａ）は実施例１における大動脈の第１ＴＤＣの一例を示した説明図であり、（Ｂ）は実施例１における心筋内冠動脈の第２ＴＤＣの一例を示した説明図である。 （Ａ）は実施例２における大動脈の第１ＴＤＣの一例を示した説明図であり、（Ｂ）は実施例２における心筋内冠動脈の第２ＴＤＣの一例を示した説明図である。
　以下、本発明が適用された実施例について図面を用いて説明する。なお、本発明の実施の形態は、下記の実施例に何ら限定されることはなく、本発明の技術的範囲に属する限り種々の形態を採りうる。
＜画像処理装置１００の構成＞
　図１を参照して、まずは、本実施例の画像処理装置１００の構成について説明する。図１は、画像処理装置１００の全体構成図である。画像処理装置１００は、例えば、汎用的なコンピュータシステムにより構成され、以下に説明する画像処理装置１００内の個々の構成要素又は機能を、コンピュータが読み取り可能な記録媒体等に保存されたコンピュータプログラムを実行することにより実現される。本実施例の画像処理装置１００は、画像データ記憶部１１０と、入力関数導出部１２０と、入力関数データ記憶部１２９と、出力関数導出部１３０と、出力関数データ記憶部１３９と、血流解析処理部１５０とを備える。また、画像処理装置１００には、入力装置１６０と、表示装置１７０と、が接続されている。
　画像処理装置１００は、被験者に造影剤を投与したあとに撮像された被験者の臓器のＣＴ画像を用いて、その臓器の血流量の定量解析を行う装置である。本実施形態では、画像処理装置１００は、被験者の心臓を心電図に同期して撮像した同位相の複数枚のＣＴのフレーム画像を対象とし、これら心臓のＣＴ画像を画素単位で解析することにより、心臓および心臓近傍の血流量の定量値を算出するために用いる入力関数及び出力関数等を特定する。本実施形態において、画像処理装置１００は、静脈に投与された造影剤によるＣＴ画像の画素値の変化に基づいて心臓等の血流量を推定し、例えば、対象となる部位に造影剤が到達することによってＣＴ画像の画素値が急激に変化した場合、画素値が変化し始める前の画素値をベース値としたり、画素値が変化し始める時刻をＡＴとしたりする。なお、本造影剤としては、種々のものを利用することができるが、本実施形態においては、例えば、イオプロミドを利用するものとする。
　画像データ記憶部１１０は、コンピュータシステムにおける一般的なメモリであって、本実施形態では、心臓の撮像によって得られた心臓及び心臓近傍のＣＴ画像の画像データを記憶する。ＣＴ画像は、造影剤が投与された後の被験者の臓器を撮像した時系列の複数フレームのＣＴ画像であって、複数のスライス画像（短軸断層画像）で構成される３次元画像（３次元ボクセルデータ）であり、詳細については後述する。
　入力関数導出部１２０は、画像データ記憶部１１０に保存されている画像データに基づいて、被験者の臓器に流入する造影剤がＣＴ画像の画素値（以下、「ＣＴ値」と称する。）に及ぼす変化に関係する入力関数等を導出する。なお、入力関数導出部１２０は、以下に説明するＲＯＩ設定部１２１と、ＲＯＩ経時変化特定部１２３と、入力関数化処理部１２５と、を有している。入力関数データ記憶部１２９は、入力関数導出部１２０で導出された入力関数等に係るデータを後述する入力関数テーブル７００として記憶する。
　ＲＯＩ設定部１２１は、医師等の解析者による入力装置１６０を介した操作に従って、ＣＴ画像において入力関数を特定するための領域であるＲＯＩ（Ｒｅｇｉｏｎ　Ｏｆ　Ｉｎｔｅｒｅｓｔ；関心領域）を設定する。具体的に、ＲＯＩ設定部１２１は、画像データ記憶部１１０から画像データを読み出して、対象とする領域が写っているスライスのフレーム画像を表示装置１７０に表示させ、解析者が入力装置３を用いて行った入力に従って、表示装置１７０に表示されている画像上にＲＯＩを設定する。なお、このＲＯＩの位置は、同一スライスの全フレームに対して共通であるものとする。
　ＲＯＩ経時変化特定部１２３は、時系列の複数フレームのＣＴ画像の画像データに基づいて、ＲＯＩ内のＣＴ値の経時変化を特定する。具体的に、ＲＯＩ経時変化特定部１２３は、ＣＴ画像における所定の領域のＣＴ値の経時変化として、ＲＯＩ内のＣＴ値から定まるＲＯＩ値のＴＤＣを生成する。なお、ＲＯＩ値とは、ＲＯＩ内のＣＴ値を統計的アルゴリズムで処理した値（統計値）であり、例えば、平均値、最頻値、中央値、最大値、最小値などのいずれかを採り得るものとする。
　入力関数化処理部１２５は、ＲＯＩ経時変化特定部１２３により生成されたＴＤＣに基づいて、ＲＯＩ内のＣＴ値が最大値Ｍａｘとなる時刻であるＭＴ、ＲＯＩが設定された領域に造影剤が到達した時刻であるＡＴ、後述する第２ＴＤＣを得るための撮像する期間の撮像区間であるＭＢ、入力関数等を特定する。なお、入力関数化処理部１２５は、ＲＯＩ経時変化特定部１２３によって生成されたＴＤＣにおいて、ＴＤＣがなめらかな曲線になるように平滑化処理を実行したうえで、入力関数等を特定する。また、入力関数化処理部１２５は、ＣＴ値が急激に上昇した後の所定のフレームである上限フレームＦａや、後述するベース値であるＢＬも特定する。
　出力関数導出部１３０は、画像データ記憶部１１０に保存されている画像データに基づいて、被験者の臓器の血管に造影剤が流入したことによりＣＴ画像のＣＴ値に及ぼす変化に関係する出力関数等を導出する。なお、出力関数導出部１３０は、ＲＯＩ単位ではなく画素単位に出力関数を導出するものであり、対象画素抽出部１３１と、ＣＴ経時変化特定部１３３と、出力関数化処理部１３７と、を有している。出力関数データ記憶部１３９は、出力関数導出部１３０で導出された出力関数に係るデータを後述する出力関数テーブル７５０として記憶する。
　対象画素抽出部１３１は、出力関数の導出対象となる画素を抽出し、フレーム画像上でその位置を特定する。具体的に、対象画素抽出部１３１は、抽出処理として、時系列の複数フレームのＣＴ画像のすべてにおいて、スライス単位にＣＴ値が所定の範囲内（例えば、５０～１６０）である画素を選択する。なお、この抽出処理において、目的とする対象画素の領域内で小さな対象外となる領域が含まれる場合には、その対象外の領域の大きさに応じて（例えば、画素数が所定数以下）、この領域を処理上の欠損とみなし、対象画素に転換する処理を実行することができる。また、対象画素の領域外で孤立した対象となる画素がある場合は、孤立した対象領域の大きさに応じて（例えば、画素数が所定数以下）、この領域を対象外とする処理を実行することもできる。そして、対象画素抽出部１３１は、このようにスライス毎に抽出された画素の位置を、出力関数を算出する対象である対象画素位置（臓器内画素位置）として特定する。
　ＣＴ経時変化特定部１３３は、時系列の複数フレームのＣＴ画像の画像データに基づいて、対象画素位置（臓器内画素位置）の画素のＣＴ値の経時変化を特定する。具体的に、ＣＴ経時変化特定部１３３は、スライス毎に、対象画素位置のＣＴ値の経時変化としてＴＤＣを生成する。
　出力関数化処理部１３７は、ＣＴ経時変化特定部１３３により生成されたＴＤＣに基づいて、出力関数等を特定する。具体的に、出力関数化処理部１３７は、ＣＴ経時変化特定部１３３で特定された経時変化であるＴＤＣに基づいて各ＣＴ値間の傾きの値を算出し、算出された傾きの値のうちで所定値以上の傾きに基づき出力関数を特定する。
　血流解析処理部１５０は、入力関数及び出力関数に基づいて、血流量の定量解析を行う。具体的に、血流解析処理部１５０は、入力関数データ記憶部１２９の入力関数テーブル７００及び出力関数データ記憶部１３９の出力関数テーブル７５０を参照して、所定の解析処理を行う。この血流量の定量解析の手法としては、例えば、Ｐａｔｌａｋ　Ｐｌｏｔ法や、Ｄｅｃｏｎｖｏｌｕｔｉｏｎ法などを用いて、ＭＢＦ（Ｍｙｏｃａｒｄｉａｌ　Ｂｌｏｏｄ　Ｆｌｏｗ：心筋血流量）、ＭＢＶ（Ｍｙｏｃａｒｄｉａｌ　Ｂｌｏｏｄ　Ｖｏｌｕｍｅ：心筋血液量）、ＭＴＴ（Ｍｅａｎ　Ｔｒａｎｓｉｔ　Ｔｉｍｅ：平均通過時間）等を特定する。具体的には、後述する入力関数テーブル７００の関数Ｆの傾きａ又は関数Ｆのデータと、出力関数テーブル７５０の傾きα又は関数Ｆのデータと、に基づいて所定の解析処理が実行される。また、血流解析処理部１５０による定量解析の結果は、表示装置１７０に表示される。このとき、例えば、定量解析の結果をＣＴ画像と並べてまたは重ねて表示して、二つの画像を比較することにより、心筋梗塞または狭心症などの際の冠動脈狭窄の評価と心筋虚血評価を同時に行うようにすることができる。また、例えば、特定されたＭＢＦと正常なＭＢＦとを比較し、特定されたＭＢＦが正常なＭＢＦよりも低い場合には、虚血状態であるものと判断でき、延いては、心筋の冠動脈における心筋梗塞または狭心症などの際の冠動脈狭窄の評価と心筋虚血評価をすることができる。
＜ＣＴ画像データの構成＞
　図２を参照して、画像データ記憶部１１０に記憶されるＣＴ画像について説明する。図２は、ＣＴ画像の一例を例示する説明図である。図２に示したＣＴ画像は、心臓を心電図に同期して撮像したものであって、心拍が同位相の所定枚数（例えば３０枚）のフレーム画像２１０である。１フレームは、複数枚のスライス画像によって構成されているとともに、１フレームのデータは、３次元の画像データで構成されている。なお、図２では、１つのスライスのフレーム画像を例示している。ＣＴ画像は、被験者に造影剤を投与した直後から撮像が開始されたものでもよく、撮像順に第１フレームから第３０フレームまでが含まれてもよい。
＜入力関数導出処理の内容＞
　図３および図４を参照して、次に、入力関数導出部１２０によって実行される入力関数導出処理について説明する。図３は、入力関数導出処理を示すフローチャートである。図４は、入力関数特定処理を示すフローチャートである。
　図３に示すように、入力関数導出処理において、まず、ＲＯＩ設定部１２０は、解析者による選択処理に沿って、画像データ記憶部１１０から特定の領域が写っているスライス画像のフレーム画像の画像データを読み出し（Ｓ３１０）、選択されたスライス画像の中で、特定の領域が鮮明に写っているフレーム画像が表示装置５に表示されているとき、解析者の操作に従ってＲＯＩを設定する（Ｓ３３０）。次に、ＲＯＩ経時変化特定部１２３は、ＲＯＩが設定された全フレーム画像についてそれぞれＲＯＩ値を算出し、ＲＯＩ値に基づくＴＤＣを描画する（Ｓ３５０）。そして、入力関数化処理部１２５は、後述する入力関数特定処理を実行する（Ｓ３９０）。
　図４に示すように、入力関数特定処理において、まず、入力関数化処理部１２５は、ＲＯＩ経時変化特定部１２３により生成されたＴＤＣに基づいて、入力関数の導出対象とするフレームの範囲を規定する上限フレームＦａを特定し（Ｓ４０５）、フレームナンバーを示す変数ｎを１に設定する（Ｓ４１０）とともに、ｎに１を加算する（Ｓ４２０）。次に、入力関数化処理部１２５は、第ｎフレームより以前のフレームのＲＯＩ値を一つの直線Ｌに近似し（Ｓ４３０）、第ｎフレーム以降で上限フレームＦａ以前のフレームのＲＯＩ値を２次関数Ｆに近似する（Ｓ４４０）。次に、入力関数化処理部１２５は、直線Ｌと第ｎフレームより以前のフレーム のＲＯＩ値の最小二乗誤差を算出するとともに、２次関数Ｆと第ｎフレーム以降で上限フレームＦａ以前のフレームのＲＯＩ値との最小二乗誤差を算出して、これらの誤差の総和を算出する（Ｓ４５０）。そして、入力関数化処理部１２５は、ステップＳ４１０～Ｓ４５０までの処理を、ｎが上限フレームＦａ－１よりも小さい間は繰り返し（Ｓ４６０でＹＥＳ）、ｎ＝２から上限フレームＦａ－１までの最小二乗誤差の総和をすべて算出する。次に、入力関数化処理部１２５は、ｎ＝２から上限フレームＦａ－１までの最小二乗誤差の総和をすべて算出すると（Ｓ４６０でＮＯ）、その中から最小二乗誤差の総和が最も小さいｎを特定する（Ｓ４７０）。そして、入力関数化処理部１２５は、ＴＤＣにおけるＲＯＩ値が最大値ＭａｘとなるフレームをＭＴとし、ステップＳ４７０で特定されたｎのときの直線Ｌ及び２次関数Ｆの交点となるフレームをＡＴとし、ＡＴとしたフレームにおけるＲＯＩ値（直線Ｌの高さ（Ｙ切片））をＢＬとし、ＭＴの前後の所定の心拍数に相当する期間を撮像区間ＭＢとするとともに、入力関数および入力関数の傾きａを特定する。
＜出力関数導出処理の内容＞
　図５および図６を参照して、出力関数導出部１３０において実行される出力関数導出処理について説明する。図５は、出力関数導出処理を示すフローチャートである。図６は、出力関数特定処理を示すフローチャートである。
　図５に示すように、出力関数導出処理において、まず、対象画素抽出部１３１は、一つのスライス画像を選択し、画像データ記憶部１１０からそのスライスのフレーム画像の画像データを読み出し（Ｓ５１０）、選択されたスライスの全画素の中から、全フレーム画像におけるＣＴ値が前述の条件を満たす画素が特定領域の画素とみなして、これらの画素を対象画素として抽出し（Ｓ５２０）、抽出された対象画素から一つの画素を選択する（Ｓ５３０）。次に、ＣＴ経時変化特定部１３３は、対象画素抽出部１３１によって選択された画素についてのＴＤＣを描画する（Ｓ５４０）。次に、出力関数化処理部１３７は、後述する出力関数特定処理をＣＴ値に基づき実行して出力関数等を特定する（Ｓ５５０）。そして、出力関数導出部１３０は、ステップＳ５３０～Ｓ５５０までの処理をすべての対象画素について行う（Ｓ５８０）とともに、ステップＳ５１０～Ｓ５８０までの処理をすべてのスライスについて行う（Ｓ５９０）。
　図６に示すように、出力関数特定処理において、まず、出力関数化処理部１３７は、ＣＴ経時変化特定部１３３によって生成されたＴＤＣにおける各プロット間の傾きαを算出し（Ｓ６１０）、算出された傾きにおいて所定値以上（例えば、０．２以上）となる傾きのみを抽出する（Ｓ６２０）。そして、出力関数化処理部１３７は、抽出された傾きαに基づき所定の関数で近似し（Ｓ６３０）、出力関数を特定する（Ｓ６５０）。
＜入力関数テーブル７００および出力関数テーブル７５０の構成＞
　図７（Ａ），（Ｂ）を参照して、入力関数データ記憶部１２９に記憶される入力関数テーブル７００および出力関数データ記憶部１３９に記憶される出力関数テーブル７５０のデータ構成について説明する。図７（Ａ），（Ｂ）は、入力関数テーブル７００及び出力関数テーブル７５０のデータ構造の一例を示した説明図であり、（Ａ）は入力関数データ記憶部１２９の入力関数テーブル７００のデータ構成の一例を示した説明図であり、（Ｂ）は出力関数データ記憶部１３９の出力関数テーブル７５０のデータ構成の一例を示した説明図である。
　図７（Ａ）に示すように、入力関数テーブル７００は、ＡＴ７１１と、ＭＴ７１２、ＢＬ７１３と、関数Ｆの傾きａ７１４と、関数Ｌ及びＦ７１５と、ＲＯＩ値７１７と、ＭＢ７１８と、をデータ項目として含んでいる。なお、関数Ｌ及びＦ７１５は、入力関数化処理部１２５で求められたｎ＝１から上限フレームＦａまでを近似した直線Ｌの式及び２次関数Ｆの式であり、ＲＯＩ値７１７は、ＴＤＣの元となるフレーム毎のＲＯＩ値である。また、血流解析処理部１５０は、入力関数テーブル７００における入力関数を血流解析に使用する場合に、ＡＴ７１１及びＢＬ７１５が原点となるようにＲＯＩ値７１７もしくは関数Ｌ及びＦ７１５を補正したものを入力関数として使用されるようにしてもよい。
　図７（Ｂ）に示すように、出力関数テーブル７５０は、スライス番号７６１と、画素位置７６２と、傾きα７６３と、関数Ｆ７６５と、をデータ項目として含んでいる。なお、画素位置７６２は、傾きα７６３を算出するための基となる２つのＣＴ値の位置を示すデータであり、関数Ｆ７６５は、傾きα７６３から算出された近似式である関数Ｆの式である。
＜心臓のＣＴ画像に対する画像処理の内容＞
　次に、本実施例における心臓のＣＴ画像の画像処理について具体的に説明する。本実施例における心臓のＣＴ画像に基づくＴＤＣに対する画像処理においては、心臓に対して、造影剤を生理食塩水で所定の濃度に希釈したもの（例えば、造影剤を生理食塩水で４倍に希釈したもの）を利用するテストインジェクション法での撮影を実施した後に、希釈しない造影剤を利用した撮影を実施する、という流れである。具体的に、まず、テストインジェクション法で撮影が実施されることによって、被験者の大動脈を含む心臓のＣＴ画像である第１ＣＴ画像が撮像される。そして、この第１ＣＴ画像に基づいて、画像処理装置１００による上述のような入力関数導出処理が実行されることによって、大動脈の第１ＴＤＣが生成されるとともに、ＭＢや入力関数等が特定される。次に、第１ＴＤＣに基づいて特定されたＭＢにおいて、希釈しない造影剤を利用した撮影が実行されることによって、被験者の大動脈および心筋内冠動脈を含む心臓のＣＴ画像である第２ＣＴ画像が撮像される。そして、この第２ＣＴ画像に基づいて、画像処理装置１００による上述のような出力関数導出処理が実行されることによって、心筋内冠動脈の第２ＴＤＣが生成されるとともに、出力関数等が特定される。そして、このように特定された入力関数や出力関数等に基づいて、画像処理装置１００によって心筋内冠動脈の血流量の定量解析が実行されることとなる。
　図８（Ａ），（Ｂ）を参照して、ここで、第１ＴＤＣおよび第２ＴＤＣに対する画像処理の概要について説明する。図８（Ａ），（Ｂ）は、本実施例におけるＴＤＣの一例を例示した説明図であり、（Ａ）は本実施例における大動脈の第１ＴＤＣの一例を示した説明図であり、（Ｂ）は本実施例における心筋内冠動脈の第２ＴＤＣの一例を示した説明図である。なお、図８（Ａ）は、心臓の大動脈近傍の第１ＣＴ画像において大動脈にＲＯＩを設定した場合のものであり、図８（Ｂ）は、対象画素を所定のＣＴ値（本実施例では、５０～１５０ＨＵとする。）に設定した場合のものである。
　図８（Ａ）に示すように、まずは、テストインジェクション法での撮影で得られた心臓の第１ＣＴ画像に基づく画像処理装置１００の画像処理によって、大動脈の第１ＴＤＣが生成されるとともに、この第１ＴＤＣに基づいてＭＢ等が特定される。具体的には、入力関数導出部１２０のＲＯＩ経時変化特定部１２３は、まず、ＲＯＩ設定部１２１で設定されたＲＯＩを、同一スライスの他のフレーム画像２１０にも適用し、フレーム画像毎に、ＲＯＩ内の画素のＣＴ値に基づいてＲＯＩを代表するＲＯＩ値を決定する。次に、ＲＯＩ経時変化特定部１２３は、同一スライスの全フレームについて定められたＲＯＩ値をプロットすることで、第１ＴＤＣを生成する。次に、ＲＯＩ経時変化特定部１２３は、第１ＴＤＣにおけるＣＴ値の最大値Ｍａｘ、またはカーブが大きく立ち上がった （ＣＴ値が急激に上昇した）後の最初のピーク値を検出する。そして、ＲＯＩ経時変化特定部１２３は、ＴＤＣにおけるＲＯＩ値の最小値Ｍｉｎから、この最大値Ｍａｘまたはピーク値までの一定割合（例えば、７０％、８０％、９０％などを設定することができ、本実施例では、７０％とする。）の値をとる。
　そして、入力関数化処理部１２５は、最大値Ｍａｘまたはピーク値が検出されたフレームよりも以前のフレームであって、上述のようにＲＯＩ経時変化特定部１２３によって特定された一定割合の値をとるフレームを上限フレームＦａとする。次に、入力関数化処理部１２５は、第ｎ（ｎは、２から上限フレームＦａ－１まで）フレームＦｎ以前のＲＯＩ値に対して最小二乗法等を適用して直線近似を行い、直線Ｌの式を導出するとともに、第ｎフレームＦｎから上限フレームＦａまでのＲＯＩ値に対して最小二乗法等を適用して２次関数に近似し、２次関数Ｆの式を導出する。次に、入力関数化処理部１２５は、直線Ｌと第ｎフレームＦｎ以前のＲＯＩ値との誤差の二乗和を算出するとともに、２次関数Ｆと第ｎフレームＦｎ以降のＲＯＩ値との誤差の二乗和（残差平方和）を算出する。そして、入力関数化処理部２５は、すべてのｎについて上記の処理を行い、これら誤差の二乗和の合計が最も小さくなるｎを特定し、この誤差が最小となるｎの場合の直線Ｌ及び２次関数Ｆをそれぞれ近似された入力関数として特定するとともに、入力関数の傾きａを特定する。ここで、ＲＯＩ値が最大値ＭａｘとなるフレームがＭＴとなり、直線Ｌ及び２次関数Ｆの交点（境界）となるフレームがＡＴとなり、直線Ｌ及び２次関数Ｆの交点（境界）となるフレームのＲＯＩ値がＢＬとなり、ＭＴの前後で３心拍に相当する区間がＭＢとなる。このように入力関数化処理部１２５は、第１ＴＤＣに基づいて、大動脈の領域のＭＴ、ＡＴ、ＢＬ、ＭＢ、傾きａ、入力関数を特定する。なお、上述のＢＬの特定において、ＲＯＩ値そのものをＢＬとする以外に、例えば、直線Ｌの高さ（Ｙ切片）をＢＬとしても良い。
　ここで、第２ＣＴ画像を得るための撮影を実施するタイミングとなる撮像区間であるＭＢについて説明する。一般的には、大動脈に造影剤が流れ込んだ後に心筋内冠動脈に造影剤が流れ込むこととなるため、ＭＴの前後であって所定の心拍数に相当する区間における撮影でＣＴ画像を取得すれば、造影剤が流れ込んだ心筋内冠動脈に対するＣＴ画像を取得することができる。特に、第２ＣＴ画像において、心筋内冠動脈に造影剤が確実に流れ込んだときのフレームにおける心筋内冠動脈のＣＴ値を得るために、本実施例における撮像区間ＭＢは、ＭＴの前の１心拍とＭＴの後の２心拍とに相当する区間としている。
　図８（Ｂ）に示すように、次に、ＭＢにおける撮影で得られた第２ＣＴ画像に基づく画像処理装置１００の画像処理によって、心筋内冠動脈の第２ＴＤＣが生成されるとともに、この第２ＴＤＣに基づいて出力関数等が特定される。ここで、心筋内冠動脈に血液が主に流れるのは、心拍動に伴う心臓の拡張期であるため、心臓の拡張期と非拡張期とではＣＴ値に変化が生じてしまい、第２ＴＤＣは段差が形成されたようなプロットとなってしまう。そこで、第２ＴＤＣに基づき出力関数を特定する場合には、所定値以上となる傾きの値に基づいて近似することによって精度の高い出力関数を特定することができる。具体的には、撮像区間ＭＢでの撮影に基づく第２ＴＤＣにおける各プロット間の傾きの値を算出して、全ての傾きの値のうち所定値以上となる傾きα１～α３の値を抽出し、この傾きα１～α３に基づく出力関数特定処理によって出力関数が特定されることとなる。
　ここで、上述のように心臓の拡張期と非拡張期との関係に沿った第２ＴＤＣにおける１心拍中の傾きの値の変化の態様は、前半は小さい傾きの値であり、後半は大きい傾きの値となるような態様である。そして、出力関数特定処理において利用する傾きの値は、１心拍中における後半の大きい傾きの値である傾きα１～α３であるものとしている。つまり、言い換えれば、出力関数特定処理は、第２ＴＤＣにおいて１心拍中の傾きのうち後半部分の傾きの値を利用した処理である。なお、出力関数を特定するときにおいて、傾きα１～α３の平均値を利用して出力関数を特定してもよいし、傾きα１～α３の各々の近似式を特定してから種々の特定方法によって出力関数を特定してもよい。
　図９（Ａ），（Ｂ）を参照して、次に、実施例２における心臓のＣＴ画像の画像処理について具体的に説明する。図９（Ａ），（Ｂ）は、実施例２におけるＴＤＣの一例を例示した説明図であり、（Ａ）は実施例２における大動脈の第１ＴＤＣの一例を示した説明図であり、（Ｂ）は実施例２における心筋内冠動脈の第２ＴＤＣの一例を示した説明図である。なお、実施例２における画像処理装置１００の構成は実施例１と同様のものであり、実施例２における入力関数導出処理および出力関数導出処理のフローは実施例１と略同様のものであるため、＜画像処理装置１００の構成＞、＜入力関数導出処理の内容＞、および＜出力関数導出処理の内容＞等の説明は省略する。
　図９に示すように、実施例２における心臓のＣＴ画像の画像処理においても、上述の実施例１と同様に、心臓に対して、テストインジェクション法での撮影を実施した後に、希釈しない造影剤を利用した撮影を実施する、という流れである。但し、上述の実施例１とは異なり、実施例２における第１ＴＤＣに基づき特定されるＭＢは１心拍に相当区間であり、第２ＴＤＣも１心拍に相当するものが生成される。そして、撮像区間ＭＢでの撮影に基づく第２ＴＤＣにおける傾きα４を抽出し、この傾きα４に基づく出力関数特定処理によって出力関数が特定されることとなる。なお、実施例２における撮像区間であるＭＢは、心筋内冠動脈に造影剤が確実に流れ込んだときのフレームにおける心筋内冠動脈のＣＴ値を得るために、ＭＴの後の１心拍に相当する区間としている。
＜実施例１および実施例２の特徴＞
　以上のように、上述の実施例１および実施例２の画像処理装置１００によれば、Ｓ５４０において、撮像区間であるＭＢでの撮影による時系列の複数フレームの第２ＣＴ画像からなる画像データに基づきＣＴ値の経時変化である第２ＴＤＣを特定し、Ｓ６２０において、第２ＴＤＣにおける各プロット間のＣＴ値の傾きであって所定の傾きα１～α３、α４を抽出し、Ｓ６３０において、これら所定の傾きα１～α３、α４に基づき所定の関数で近似することができる。特に、撮像区間であるＭＢにおける算出された複数の傾きに対して、所定値より大きい傾きα１～α３、α４を抽出した処理を実行することができる。したがって、このような実施例１および実施例２の画像処理装置１００であれば、例えば、心筋内冠動脈のＣＴ画像に対する画像処理に利用すると、本来、心筋血流が流れていない心臓の非拡張期のデータを排除して、拡張期のデータに基づいた精度の高い画像処理、延いては、精度の高い血流量の定量解析を実行することができる。また、従来の画像処理に比べて簡便な方法で解析することができるため、画像処理による画像処理装置１００における負荷を軽減することもできる。
　また、上述の実施例１および実施例２の画像処理装置１００によれば、画像データ記憶部１１０は、第１所定数のフレームの第１ＣＴ画像からなる第１画像データと、第１所定数に比べて少ない第２所定数のフレームの第２ＣＴ画像からなる第２画像データと、を記憶しており、希釈した造影剤を利用したテストインジェクション法の撮影による第１ＣＴ画像から特定されたＭＴに基づいて、所定の心拍数に相当する撮像区間であるＭＢを特定し、ＭＢでの希釈していない造影剤を利用した撮影が実施されて、得られた第２ＣＴ画像に基づいて、画像処理装置１００による画像処理を実行することができる。したがって、従来のように２０～３０心拍に応じた期間のＣＴＰ検査をする場合に比べて、少ない心拍数に応じた期間である撮像区間であるＭＢでの撮影によるＣＴ画像で画像処理を実行することができるため、従来に比べて被験者における被ばく量を軽減するとともに簡便な検査に基づく画像処理を実行することができる。
　さらに、上述の実施例２の画像処理装置１００によれば、撮像区間であるＭＢに基づいた１回の拍動に相当するフレーム数の第２ＣＴ画像および第２ＴＤＣに基づいて、画像処理装置１００による画像処理を実行することができる。したがって、従来のように２０～３０心拍に応じた期間のＣＴＰ検査をする場合に比べて、少ない心拍数に応じた期間である撮像区間ＭＢのＣＴＰ検査によるＣＴ画像で画像処理を実行することができるため、従来に比べて被験者における被ばく量を軽減するとともに簡便な検査に基づく画像処理を実行することができる。
＜その他の実施例＞
　上述の実施例において、第１ＣＴ画像に基づくＭＴの前後で３心拍又は１心拍に相当する区間を撮像区間であるＭＢとしたが、本発明ではこれに限定されず、ＭＢを以下のようにしてもよい。例えば、ＭＴの前後で２心拍や４～１０心拍に相当する期間を撮像区間であるＭＢとしてもよい。このようなＭＢや撮像タイミングに基づく撮影であっても、従来に比べて被験者における被ばく量を軽減するとともに簡便な検査に基づく画像処理を実行することができる。
　上述の実施例において、図６に示した出力関数特定処理では、算出された傾きにおいて０．２以上となる傾きのみを抽出するものとしたが、本発明ではこれに限定されず、０．２以外の数値を閾値としてもよい。例えば、複数回の撮影による結果から算出された統計値に基づく数値を閾値としてもよい。このような閾値であれば、より精度の高い画像処理、延いては、精度の高い血流量の定量解析を実行することができる。
　上述の実施例において、図６に示した出力関数特定処理では、算出された傾きにおいて０．２以上となる傾きのみを抽出し、つまり、心臓の拡張期における傾き（１心拍中の傾きのうち後半の傾き）の値のみを抽出するものとしたが、本発明ではこれに限定されず、心臓の非拡張期におけるＣＴ値の変化に基づき出力関数を特定するようにしてもよい。具体的には、例えば、第２ＴＤＣにおける０．２未満となる傾きを抽出し、抽出された傾きに基づいた各ＣＴ値の変化（心臓の非拡張期におけるＣＴ値の変化）の差分に基づく差分法等によって出力関数を特定する、といった処理であってもよい。
　上述の実施例において、図６に示した出力関数特定処理は、出力関数導出処理において利用するようにしたが、本発明はこれに限定されず、入力関数導出処理において図６に示した出力関数特定処理のような処理を利用してもよい。例えば、第１ＴＤＣに基づく入力関数を特定するときに出力関数特定処理のように傾きから入力関数を特定するようにしてもよい。このような画像処理であれば、従来の画像処理に比べてより簡便な方法で解析することができるため、画像処理による画像処理装置１００における負荷をより軽減することができる。
　上述の実施例において、第１ＣＴ画像を撮像するために、造影剤を生理食塩水で希釈した造影剤を利用したテストインジェクション法でＤｙｎａｍｉｃ撮影を実施していたが、本発明ではこれに限定されない。例えば、対象部位にＲＯＩを設定してリアルタイムでＲＯＩのＣＴ値（ＲＯＩ値）をモニタリングするボーラストラッキング法で撮影を実施してもよい。また、希釈していない造影剤を利用したテストインジェクション法でＤｙｎａｍｉｃ撮影を実施してもよく、さらに、テストインジェクション法とボーラストラッキング法とを組み合わせた、所謂、テストボーラストラッキング法でＤｙｎａｍｉｃ撮影を実施してもよい。つまり、上述の実施例において、第２ＣＴ画像を撮像するための撮影を実施するタイミングを最適化するためには、種々の手法を利用して第１ＣＴ画像を撮像することができる。
　上述の実施例においては、希釈した造影剤を利用したテストインジェクション法での撮影で第１ＣＴ画像を得て、その後に希釈していない造影剤を利用した撮影で第２ＣＴ画像を得る、といった手順であるものとしたが、本発明はこれに限定されない。例えば、第１ＣＴ画像を撮像するためのテストインジェクション法での撮影の実行中に、所定の心拍においてのみ複数心位相を撮影することによって第２ＣＴ画像を撮像する、といった手法であってもよい。
　上述の実施例において、対象画素抽出部１３１は、上述のような抽出処理を実行していたが、本発明ではこれに限定されず、画素単位に以下の処理を実行してもよい。例えば、対象画素抽出部１３１が、時系列の複数フレームのＣＴ画像において、各画素のＣＴ値の変化量、例えば、最大値と最小値との差分の値に基づいて画素を選択してもよい。具体的に、対象画素抽出部１３１は、一つのスライスにおいて、全フレーム画像の中でＣＴ値の最大値と最小値との差を求め、その差が所定の値（例えば５０～１５０）である画素を抽出してもよい。さらに、対象画素抽出部１３１は、このような抽出処理を上述の抽出処理に加えて実行してもよいし、いずれか一方のみを実行してもよい。
　上述の実施例において、入力関数特定処理では上限フレームＦａを上述のような方法で特定していたが、本発明はこれに限定されない。例えば、第１ＴＤＣの変化率を求めて、変化率が最大となるフレームと変化率が０となるフレームとの間のいずれかのフレームを上限フレームＦａとしてもよい。
　上述の実施例において、入力関数特定処理では直線と２次曲線との２つの関数で第１ＴＤＣを関数化しているが、本発明はこれに限定されない。例えば、第ｎフレームＦｎ以前を直線Ｌで近似し、第ｎフレームＦｎ以降のＲＯＩ値又はＣＴ値の分布を直線または３次以上の多次関数で近似してもよい。あるいは、上限フレームＦａまでのすべてのＲＯＩ値又はＣＴ値の分布を多次元多項式で表される関数に近似してもよい。また、第１ＴＤＣを３つ以上の関数で近似してもよく、例えば、第２フレームから上限フレームＦａ－１までを３つ以上の区間に分けて、それぞれの区間を所定の関数で近似してもよい。第１ＴＤＣにおいては、カーブが大きく立ち上がる（ＲＯＩ値又はＣＴ値が急激に上昇する）直前に、一旦ＲＯＩ値又はＣＴ値が低下する場合がある。このような現象が現れた場合には、ＲＯＩ値又はＣＴ値がほぼ一定で直線近似できる区間（第１区間）と、ＲＯＩ値又はＣＴ値が減少する区間（第２区間）と、その後のＲＯＩ値又はＣＴ値が急激に上昇する区間（第３区間）とに分けてもよい。そして、例えば、第１区間を直線、第２区間を２次以上の関数、第３区間を別の２次以上の関数で近似するようにしてもよい。このような場合には、第１区間と第２区間の境界となるフレーム、または第２区間と第３区間の境界となるフレームの何れかをＡＴとし、このＡＴとなるフレームのＲＯＩ値又はＣＴ値をベース値としてもよい。
　さらに、第１ＴＤＣを一つの関数で近似するようにしてもよく、例えば、正規累積分布関数または累積分布関数を第１ＴＤＣにフィッティングさせて、第１ＴＤＣを単一の関数で近似してもよい。このように第１ＴＤＣを正規累積分布関数で関数化する場合、第１ＴＤＣの立ち上がりのカーブに最もよくフィットするように正規分布の標準偏差（ＳＤ）及び平均値を選択するようにしてもよい。このとき、例えば、第１ＴＤＣに近似された正規累積分布関数の－３ＳＤに最も近いフレームをＡＴとし、このＡＴとなるフレームのＲＯＩ値又はＣＴ値をベース値としてもよい。
　上述の実施例において、関数特定処理では平滑化した第１ＴＤＣに対して関数近似していたが、本発明はこれに限定されず、上述のようにｍ次関数にフィッティングする以外の方法で第１ＴＤＣが平滑化されてもよい。例えば、第１ＴＤＣを生成するときに、周辺の画素との平均値を求めてこれを利用して平滑化してもよいし、第１ＴＤＣを移動平均によって平滑化してもよい。また、平滑化していない第１ＴＤＣに対して関数近似してもよい。
　上述の実施例において、上述のように血流解析処理部１５０による定量解析の結果が表示装置１７０に表示されていたが、本発明はこれに限定されない。例えば、心臓領域を任意の複数のセグメントに分けて表示する画像を表示装置１７０に表示させるとき、各セグメントはそれぞれの血流量に応じた表示態様で表示されるようにしてもよい。あるいは、画像データ記憶部１１０に保存されている画像データに基づいた臓器の３Ｄ画像を表示させるときは、その３Ｄ画像における各画素の表示態様が血流量に応じたものになるようにしてもよい。なお、ＣＴ画像は、個々の画素の座標情報を基に３Ｄ画像表示を行ってもよい。
　以上、実施例、変形例に基づき本発明について説明してきたが、上記した発明の実施の形態は、本発明の理解を容易にするためのものであり、本発明を限定するものではない。本発明は、その趣旨並びに特許請求の範囲を逸脱することなく、変更、改良され得ると共に、本発明にはその等価物が含まれる。
　１００…画像処理装置、１１０…画像データ記憶部、１２０…入力関数導出部、１２１…ＲＯＩ設定部、１２３…経時変化特定部、１２５…関数化処理部、１２９…入力関数データ記憶部、１３０…出力関数導出部、１３１…対象画素抽出部、１３３…ＣＴ経時変化特定部、１３７…出力関数化処理部、１３９…出力関数データ記憶部、１５０…血流解析処理部、１６０…入力装置、１７０…表示装置。

Claims (6)

  1. 　造影剤が投与された後の被験者の臓器を撮像したものであって、時系列に沿った所定数のフレームのＣＴ（Ｃｏｍｐｕｔｅｄ　Ｔｏｍｏｇｒａｐｈｙ）画像からなる画像データを記憶する記憶手段を備えた画像処理装置のためのコンピュータプログラムにおいて、
    　前記画像処理装置に、
    　　前記所定数のフレームのＣＴ画像からなる画像データに基づいて、ＣＴ値の経時変化を特定する第１ステップと、
    　　前記第１ステップで特定された前記ＣＴ値の経時変化に基づいて、所定時間に対する前記ＣＴ値の傾きであって所定の傾きを特定する第２ステップと、
    　　前記第２ステップで特定された前記所定の傾きに基づいて、前記ＣＴ値の経時変化を所定の関数で近似する第３ステップと、
    　を実行させる
    　ことを特徴とするコンピュータプログラム。
  2. 　前記所定の傾きは、前記所定時間に対する前記ＣＴ値の傾きであって所定値より大きい傾きである
    　ことを特徴とする請求項１に記載のコンピュータプログラム。
  3. 　前記記憶手段は、
    　　第１所定数のフレームのＣＴ画像からなる第１画像データと、
    　　前記第１所定数に比べて少ない第２所定数のフレームのＣＴ画像からなる第２画像データと、
    　を記憶しており、
    　前記画像処理装置に、
    　　前記第１画像データに基づき前記第１ステップで特定された前記ＣＴ値の経時変化において、前記ＣＴ値が最大値となるタイミングである最大タイミングを特定する第４ステップ、
    　をさらに実行させ、
    　　前記第４ステップで特定された前記最大タイミングに基づく前記第２所定数のフレームのＣＴ画像からなる前記第２画像データに基づいて、前記第１ステップ、前記第２ステップ、および前記第３ステップを実行させる
    　ことを特徴とする請求項１又は請求項２に記載のコンピュータプログラム。
  4. 　前記画像処理装置に、
    　　前記第４ステップで特定された前記最大タイミングに基づいた前記臓器の１回の拍動に相当する前記第２所定数のフレームのＣＴ画像からなる前記第２画像データに基づいて、前記第１ステップ、前記第２ステップ、および前記第３ステップを実行させる
    　ことを特徴とする請求項３に記載のコンピュータプログラム。
  5. 　造影剤が投与された後の被験者の臓器を撮像したものであって、時系列に沿った所定数のフレームのＣＴ（Ｃｏｍｐｕｔｅｄ　Ｔｏｍｏｇｒａｐｈｙ）画像からなる画像データを記憶する記憶手段と、
    　前記所定数のフレームのＣＴ画像からなる画像データに基づいて、ＣＴ値の経時変化を特定する第１手段と、
    　前記第１手段によって特定された前記ＣＴ値の経時変化に基づいて、所定時間に対する前記ＣＴ値の傾きであって所定の傾きを特定する第２手段と、
    　前記第２手段によって特定された前記所定の傾きに基づいて、前記ＣＴ値の経時変化を所定の関数で近似する第３手段と、
    　を備えることを特徴とする画像処理装置。
  6. 　造影剤が投与された後の被験者の臓器を撮像したものであって、時系列に沿った所定数のフレームのＣＴ（Ｃｏｍｐｕｔｅｄ　Ｔｏｍｏｇｒａｐｈｙ）画像からなる画像データを記憶する記憶手段を備えた画像処理装置が実行する画像処理方法であって、
    　前記画像処理装置が、
    　　前記所定数のフレームのＣＴ画像からなる画像データに基づいて、ＣＴ値の経時変化を特定する第１ステップと、
    　　前記第１ステップで特定された前記ＣＴ値の経時変化に基づいて、所定時間に対する前記ＣＴ値の傾きであって所定の傾きを特定する第２ステップと、
    　　前記第２ステップで特定された前記所定の傾きに基づいて、前記ＣＴ値の経時変化を所定の関数で近似する第３ステップと、
    　を実行することを特徴とする画像処理方法。
PCT/JP2019/021147 2018-05-29 2019-05-28 コンピュータプログラム、画像処理装置、および画像処理方法 Ceased WO2019230738A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
CN201980035304.2A CN112165901B (zh) 2018-05-29 2019-05-28 计算机可读记录介质、图像处理装置、以及图像处理方法
EP19812179.0A EP3811869B1 (en) 2018-05-29 2019-05-28 Computer program, image processing device, and image processing method
JP2020522222A JP7236747B2 (ja) 2018-05-29 2019-05-28 コンピュータプログラム、画像処理装置、および画像処理方法
US17/058,694 US11813106B2 (en) 2018-05-29 2019-05-28 Image processing device, and image processing method utilizing time-series computed tomography (CT) images

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2018102516 2018-05-29
JP2018-102516 2018-05-29

Publications (1)

Publication Number Publication Date
WO2019230738A1 true WO2019230738A1 (ja) 2019-12-05

Family

ID=68698858

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2019/021147 Ceased WO2019230738A1 (ja) 2018-05-29 2019-05-28 コンピュータプログラム、画像処理装置、および画像処理方法

Country Status (5)

Country Link
US (1) US11813106B2 (ja)
EP (1) EP3811869B1 (ja)
JP (1) JP7236747B2 (ja)
CN (1) CN112165901B (ja)
WO (1) WO2019230738A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112053413A (zh) * 2020-08-31 2020-12-08 杭州脉流科技有限公司 脑部ct灌注成像的参数图获取方法、装置、计算机设备和存储介质
JP2023539301A (ja) * 2020-08-26 2023-09-13 ロンドン ヘルス サイエンシーズ センター リサーチ インコーポレイテッド 血流イメージング

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013057982A1 (ja) * 2011-10-19 2013-04-25 株式会社日立メディコ 画像診断装置、および画像判別方法
US20140163403A1 (en) * 2012-12-12 2014-06-12 The Texas A&M University System Automated determination of arterial input function areas in perfusion analysis
WO2016009957A1 (ja) 2014-07-15 2016-01-21 富士フイルムＲｉファーマ株式会社 コンピュータプログラム、画像処理装置及び方法

Family Cites Families (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7022077B2 (en) * 2000-11-28 2006-04-04 Allez Physionix Ltd. Systems and methods for making noninvasive assessments of cardiac tissue and parameters
WO2004002305A2 (en) * 2002-07-01 2004-01-08 Allez Physionix Systems and methods for making noninvasive assessments of cardiac tissue and parameters
JP5536974B2 (ja) * 2006-11-08 2014-07-02 株式会社東芝 Ｘ線診断装置及び画像処理装置
JP5010375B2 (ja) * 2007-07-18 2012-08-29 株式会社東芝 医用画像診断装置
US20100016731A1 (en) * 2008-07-15 2010-01-21 Cardiox Corporation Hemodynamic Detection of Circulatory Anomalies
CN101658428A (zh) * 2008-08-28 2010-03-03 西门子（中国）有限公司 灌注成像技术中对位图进行处理的方法和系统
US20110313287A1 (en) 2009-03-04 2011-12-22 Nemoto Kyorindo Co., Ltd. Chemical liquid injector and x-ray ct system
US8538106B2 (en) * 2009-10-20 2013-09-17 Siemens Aktiengesellschaft Three-dimensional esophageal reconstruction
JP5422543B2 (ja) * 2010-12-15 2014-02-19 富士フイルムＲｉファーマ株式会社 血流量推定装置、方法及びコンピュータプログラム
JP5753697B2 (ja) * 2011-02-07 2015-07-22 株式会社日立メディコ Ｘ線ｃｔ装置
JP6298811B2 (ja) 2013-04-11 2018-03-20 株式会社根本杏林堂 薬液注入装置
JP5695140B2 (ja) * 2013-07-29 2015-04-01 株式会社東芝 医用画像診断装置
JP6285304B2 (ja) * 2014-07-15 2018-02-28 富士フイルムＲｉファーマ株式会社 コンピュータプログラム、画像処理装置及び方法
US9517042B2 (en) * 2014-09-12 2016-12-13 General Electric Company Systems and methods for imaging phase selection for computed tomography imaging
KR102389347B1 (ko) * 2015-02-05 2022-04-22 삼성메디슨 주식회사 초음파 진단장치 및 그에 따른 초음파 진단 장치의 동작 방법
US10687700B2 (en) * 2015-08-03 2020-06-23 National University Corporation Ehime University Visual field measuring method, visual field measuring apparatus, and optotype
JP2017074123A (ja) * 2015-10-13 2017-04-20 東芝メディカルシステムズ株式会社 医用画像処理装置及びｘ線診断装置
KR101798083B1 (ko) * 2016-02-04 2017-12-12 삼성전자주식회사 단층 영상 처리 장치, 방법, 및 그 방법에 관련된 기록매체
JP6846758B2 (ja) * 2017-03-09 2021-03-24 キヤノンメディカルシステムズ株式会社 画像処理装置
JP6945462B2 (ja) * 2018-01-17 2021-10-06 富士フイルム株式会社 画像処理装置、画像処理方法、及び画像処理プログラム
JP7033931B2 (ja) * 2018-01-17 2022-03-11 富士フイルム株式会社 画像処理装置、画像処理方法、及び画像処理プログラム

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013057982A1 (ja) * 2011-10-19 2013-04-25 株式会社日立メディコ 画像診断装置、および画像判別方法
US20140163403A1 (en) * 2012-12-12 2014-06-12 The Texas A&M University System Automated determination of arterial input function areas in perfusion analysis
WO2016009957A1 (ja) 2014-07-15 2016-01-21 富士フイルムＲｉファーマ株式会社 コンピュータプログラム、画像処理装置及び方法

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2023539301A (ja) * 2020-08-26 2023-09-13 ロンドン ヘルス サイエンシーズ センター リサーチ インコーポレイテッド 血流イメージング
CN112053413A (zh) * 2020-08-31 2020-12-08 杭州脉流科技有限公司 脑部ct灌注成像的参数图获取方法、装置、计算机设备和存储介质
CN112053413B (zh) * 2020-08-31 2023-08-29 杭州脉流科技有限公司 脑部ct灌注成像的参数图获取方法、装置、计算机设备和存储介质

Also Published As

Publication number Publication date
CN112165901B (zh) 2024-04-05
EP3811869A4 (en) 2022-04-27
EP3811869B1 (en) 2024-07-03
US11813106B2 (en) 2023-11-14
JP7236747B2 (ja) 2023-03-10
JPWO2019230738A1 (ja) 2021-07-29
US20210196220A1 (en) 2021-07-01
EP3811869A1 (en) 2021-04-28
CN112165901A (zh) 2021-01-01

Similar Documents

Publication Publication Date Title
JP5851553B2 (ja) イメージング法における対象領域の識別と時間値曲線の抽出
JP5558672B2 (ja) 画像処理装置及びｘ線コンピュータ断層撮影装置
US9275266B2 (en) Apparatus and method for tracking contour of moving object, and apparatus and method for analyzing myocardial motion
JP5093727B2 (ja) 連続ｘ線画像スクリーニング検査装置、プログラム及び記録媒体
JP6196309B2 (ja) 弁を通る動くオブジェクトへの逆流フローを評価するデータを提供する方法と装置
JP6667999B2 (ja) 画像処理装置、画像処理方法、及びプログラム
KR101468418B1 (ko) 초음파 영상 처리 방법 및 장치
CN105022719B (zh) 医学造影系统及方法
US10499871B2 (en) X-ray diagnostic apparatus and image processing apparatus
US10438364B2 (en) System and method for registration of perfusion images
CN110312475A (zh) 用于对人心脏和心房进行建模的方法和系统
JP7236747B2 (ja) コンピュータプログラム、画像処理装置、および画像処理方法
JP2011043489A (ja) 心筋の局所壁運動状態を評価するための、画像処理装置、画像処理プログラム及び画像処理方法
US20110026775A1 (en) Method for measurement of a flow in an object, especially a lumen or a vessel
JP6343004B2 (ja) コンピュータプログラム、画像処理装置及び方法
JP2015217113A (ja) 血管解析装置、医用画像診断装置、血管解析方法及び血管解析プログラム
JP2003299646A (ja) 画像解析装置
CN105007818B (zh) X射线诊断装置以及图像处理装置
CN117412712A (zh) 对左心室顺应性的无创测量
JP6675458B2 (ja) 血管解析装置、血管解析方法及び血管解析プログラム
KR102000615B1 (ko) 관상동맥의 시작점을 자동으로 추출하기 위한 방법 및 그 장치
JP2014196990A (ja) Ｒｏｉ設定装置、ｒｏｉ設定方法、コンピュータプログラム及び脳血流推定装置
Zwirn et al. Automatic endocardial-boundary detection in low mechanical-index contrast echocardiography
Mor-Avi et al. Segmental analysis of Color Kinesis: new echocardiographic technique for objective assessment of regional LV endocardial motion

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: 19812179

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2020522222

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2019812179

Country of ref document: EP

Effective date: 20201201