WO2018189815A1 - 光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置 - Google Patents

光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置 Download PDF

Info

Publication number
WO2018189815A1
WO2018189815A1 PCT/JP2017/014869 JP2017014869W WO2018189815A1 WO 2018189815 A1 WO2018189815 A1 WO 2018189815A1 JP 2017014869 W JP2017014869 W JP 2017014869W WO 2018189815 A1 WO2018189815 A1 WO 2018189815A1
Authority
WO
WIPO (PCT)
Prior art keywords
property value
optical property
value distribution
measured
light
Prior art date
Application number
PCT/JP2017/014869
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 PCT/JP2017/014869 priority Critical patent/WO2018189815A1/ja
Publication of WO2018189815A1 publication Critical patent/WO2018189815A1/ja
Priority to US16/594,166 priority patent/US10955337B2/en

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/47—Scattering, i.e. diffuse reflection
    • G01N21/4795—Scattering, i.e. diffuse reflection spatially resolved investigating of object in scattering medium
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/47—Scattering, i.e. diffuse reflection
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B10/00—Other methods or instruments for diagnosis, e.g. instruments for taking a cell sample, for biopsy, for vaccination diagnosis; Sex determination; Ovulation-period determination; Throat striking implements
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/47—Scattering, i.e. diffuse reflection
    • G01N2021/4704—Angular selective
    • G01N2021/4709—Backscatter
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2201/00—Features of devices classified in G01N21/00
    • G01N2201/06—Illumination; Optics
    • G01N2201/061—Sources
    • G01N2201/06113—Coherent sources; lasers
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2201/00—Features of devices classified in G01N21/00
    • G01N2201/12—Circuits of general importance; Signal processing
    • G01N2201/126—Microprocessor processing

Definitions

  • the present invention relates to an optical property value distribution estimation method, an optical property value distribution estimation program, and an optical property value distribution estimation apparatus.
  • a method using a reflection type diffuse optical tomography method is known as a non-invasive method for estimating the distribution of optical property values inside a living tissue (see, for example, Patent Documents 1 to 3).
  • this method first, light is irradiated on the object to be measured, backscattered light from the object to be measured is received, and a measurement value of the backscattered light is acquired.
  • a numerical calculation is performed using a light propagation model assuming an optical property value distribution inside the measurement object, and backscattered light from the measurement object is calculated as a predicted value.
  • the measured object model is divided into meshes (or voxels), the scattering coefficient and the absorption coefficient are set as parameters for each mesh, and the state of light propagating inside the measured object is expressed by the light transport equation. And the light diffusion equation. Then, the measured value and the predicted value are compared, and it is determined whether or not the degree of coincidence is a predetermined value or more.
  • the optical property value distribution inside the object to be measured is assumed again, a predicted value is calculated, and this calculation is repeated until the degree of coincidence becomes equal to or higher than the predetermined value.
  • This iterative calculation is also called an inverse analysis operation.
  • the degree of coincidence can be evaluated by setting an error function from the difference between the measured value and the predicted value.
  • the difference between the measured value and the predicted value can be reduced by assuming again the optical property value distribution inside the object to be measured so that the error function becomes small.
  • Various known optimization algorithms can be applied to the inverse analysis operation described above.
  • Patent Document 1 illumination is obliquely applied to a measurement object, and backscattered light from the measurement object is received in an oblique direction opposite to the illumination direction. It is disclosed that a physical property value distribution can be estimated with high accuracy.
  • Patent Document 2 discloses that the optical property value distribution of the shallow layer of the measurement object can be estimated with higher accuracy by performing measurement in a plurality of illumination directions and a plurality of light receiving directions.
  • Patent Document 3 discloses that the backscattered light of the light reaching the deep layer of the living body is isotropic because of the large number of scattering. In other words, the backscattered light of the light reaching the deep layer of the living body has a low angle dependency of the light intensity. On the other hand, the phase shift, which is TOF (Time Of Flight) information of backscattered light from the object to be measured, tends to increase as the light reaches the deep layer. It is disclosed that an optical property value distribution can be estimated with high accuracy.
  • TOF Time Of Flight
  • optical property value distribution distributed of scattering coefficient, absorption coefficient, etc.
  • Patent Documents 1 to 3 disclose measurement methods suitable for estimation of the shallow and deep layers of biological tissue, the optical physical property value distribution can be obtained with high accuracy over a wide range from the shallow layer to the deep layer. There was a problem that it could not be estimated.
  • the present invention has been made in view of the above, and an optical property value distribution estimation method capable of estimating an optical property value distribution with high accuracy over a wide range from a shallow layer to a deep layer of biological tissue.
  • An object of the present invention is to provide an optical property value distribution estimation program and an optical property value distribution estimation device.
  • an optical property value distribution estimation method includes an optical property value distribution estimation device that estimates an optical property value distribution inside a measurement object.
  • a method for estimating an optical property value distribution to be performed wherein a first measurement value obtained by measuring isotropic backscattered light of light irradiated on the object to be measured is read from a storage unit, and the object to be measured is subjected to inverse analysis calculation.
  • the second measured value that is measured is read from the storage unit, and the second physical property value distribution inside the measured object is obtained by inverse analysis calculation using at least a part of the first optical property value distribution as an initial value.
  • the second optical property value distribution is a shallower optical property value distribution than the first optical property value inside the object to be measured. It is characterized by that.
  • the object to be measured along the light irradiated on the object to be measured in the second estimation step, in the first optical property value distribution, the object to be measured along the light irradiated on the object to be measured.
  • the second optical property value distribution is estimated by using, as an initial value, a value obtained by replacing the optical property value of the entire region of the object to be measured with the optical property value of the vertical region in the direction perpendicular to the surface of the object to be measured.
  • the first measured value is measured by irradiating the object to be measured with pulsed light or light whose intensity is periodically modulated. It is a value including information, and the first estimating step estimates the first optical property value distribution based on the TOF information.
  • the second measurement value is an angle at which light is applied to the object to be measured, or backscattered light from the object to be measured is received. It is a value obtained by measuring at least one of the angles at an acute angle with respect to the surface of the object to be measured, and the second estimation step is based on the light intensity angle distribution of the second measurement value.
  • the second optical property value distribution is estimated.
  • the first measured value is a value measured using at least one of a photoacoustic phenomenon and an orthogonal Nicol method.
  • the second measured value is a value measured using at least one of the OCT method and the parallel Nicol method. To do.
  • the second measurement value is a value measured by a measurement method determined according to an optical property value representing a known scattering intensity. It is characterized by being.
  • the measurement method predicts a depth of a region having a small optical property value in a surface layer region of the object to be measured, and according to the depth. Then, an angle for irradiating light to the object to be measured and an angle for receiving backscattered light from the object to be measured are determined, and the second measurement value is measured at the determined angle.
  • the depth is calculated by 1 / ⁇ s ′ round using an approximate equivalent scattering coefficient ⁇ s ′ round of the object to be measured. It is characterized by.
  • An optical property value distribution estimation method is an optical property value distribution estimation method performed by an optical property value distribution estimation device that estimates an optical property value distribution inside a measurement object.
  • a first measurement value obtained by measuring backscattered light of the light irradiated on the object to be measured is read from the storage unit, and includes a region that is estimated to have a large loss of light energy inside the object to be measured by inverse analysis calculation.
  • Second light which is optical property value distribution of region Characterized in that it comprises a second estimation step of estimating a physical property value distribution, the.
  • An optical property value distribution estimation method is an optical property value distribution estimation method performed by an optical property value distribution estimation device that estimates an optical property value distribution inside a measurement object.
  • the first measurement value obtained by measuring the backscattered light of the light irradiated on the first region of the object to be measured is read from the storage unit, and is the first optical property value distribution inside the object to be measured by inverse analysis calculation.
  • a second optical physical property distribution which is a shallower optical physical property value distribution than the first optical physical property value distribution inside the object to be measured by inverse analysis calculation using at least a part of the first optical physical property value distribution as an initial value.
  • An optical property value distribution estimation program is an optical property value distribution estimation program performed by an optical property value distribution estimation device that estimates an optical property value distribution inside a measurement object, A first measurement value obtained by measuring isotropic backscattered light of the light irradiated on the object to be measured is read from the storage unit, and the first optical value is an optical property value distribution inside the object to be measured by inverse analysis calculation. A first estimation step for estimating a physical property value distribution, and a second measurement value obtained by measuring backscattered light that is more anisotropic than the first measurement value of the light irradiated on the object to be measured are read from the storage unit.
  • the second optical property value distribution is a shallower optical property value distribution than the first optical property value distribution inside the object to be measured by inverse analysis calculation using at least a part of the first optical property value distribution as an initial value.
  • Second estimation step for estimating optical property value distribution And characterized by causing a computer to execute the.
  • An optical property value distribution estimation apparatus is an optical property value distribution estimation apparatus that estimates an optical property value distribution inside a measurement object, such as light irradiated on the measurement object.
  • a first measurement value obtained by measuring the anisotropic backscattered light is read from the storage unit, a first optical property value distribution which is an optical property value distribution inside the measurement object is estimated by inverse analysis calculation, and the measurement object
  • a second measurement value obtained by measuring backscattered light that is more anisotropic than the first measurement value of the light irradiated to the object is read from the storage unit, and at least a part of the first optical property value distribution is an initial value.
  • a reverse analysis calculation unit for estimating a second optical property value distribution which is a shallow optical property value distribution from the first optical property value distribution inside the object to be measured.
  • an optical property value distribution estimation method capable of accurately estimating an optical property value distribution over a wide range from a shallow layer to a deep layer of a biological tissue, and An optical property value distribution estimation device can be realized.
  • FIG. 1 is a block diagram showing a configuration of an optical property value distribution estimation apparatus according to Embodiment 1.
  • FIG. 2 is a diagram illustrating a state in which the TOF sensor unit, the angle distribution measurement illumination unit, and the angle distribution measurement light receiving unit illustrated in FIG. 1 are arranged with respect to an object to be measured.
  • FIG. 3 is a diagram illustrating a state in which the first measurement value is measured by the TOF sensor unit illustrated in FIG.
  • FIG. 4 is a diagram illustrating a state in which the second measurement value is measured by the angle distribution measuring apparatus illustrated in FIG.
  • FIG. 5 is a flowchart showing an operation of measuring an object to be measured by the optical property value distribution estimating apparatus shown in FIG.
  • FIG. 6 is a diagram illustrating light propagation in the deep layer and the shallow layer.
  • FIG. 7 is a diagram illustrating an object to be measured that is measured by the optical property value distribution estimating apparatus according to the embodiment.
  • FIG. 8 is a diagram showing the optical property values of the measurement object shown in FIG.
  • FIG. 9 is a diagram illustrating a measurement result by the TOF measurement unit.
  • FIG. 10 is a diagram illustrating a state in which the object to be measured illustrated in FIG. 7 is divided into meshes.
  • FIG. 11 is an enlarged view of the mesh shown in FIG.
  • FIG. 12 is a diagram illustrating a state in which the illumination unit and the light receiving unit are arranged on the measurement object illustrated in FIG. 10.
  • FIG. 13 is a diagram illustrating an estimation result by the first estimation step.
  • FIG. 14 is a diagram for explaining an initial value setting method according to the second estimation step.
  • FIG. 14 is a diagram for explaining an initial value setting method according to the second estimation step.
  • FIG. 15 is a diagram illustrating a measurement result by the angle distribution measurement unit.
  • FIG. 16 is a diagram illustrating a measurement result by the angle distribution measurement unit.
  • FIG. 17 is a diagram illustrating a measurement result by the angle distribution measurement unit.
  • FIG. 18 is a diagram illustrating a measurement result by the angle distribution measurement unit.
  • FIG. 19 is a diagram illustrating an estimation result by the second estimation step.
  • FIG. 20 is a diagram illustrating an example of a method for estimating the thickness of the first layer.
  • an optical property value distribution estimation method for estimating an optical property value distribution using TOF information and a light intensity angle distribution as a measurement value will be described as an example.
  • the present invention can also be applied to an optical property value distribution estimation method using values.
  • FIG. 1 is a block diagram showing a configuration of an optical property value distribution estimation apparatus according to Embodiment 1.
  • the optical property value distribution estimation apparatus 1 according to the first embodiment is an apparatus that estimates an optical property value distribution inside a measurement object. As shown in FIG. 1, the optical property value distribution estimation apparatus 1 irradiates the object to be measured with light whose intensity is periodically modulated and receives backscattered light from the object to be measured. And an angle distribution measuring apparatus 20 that irradiates light to the object to be measured with a variable angle and receives backscattered light from the object to be measured with a variable angle, and the TOF sensor unit 10 and the angle distribution measuring apparatus 20. Includes a processing device 30 that processes the signal acquired by the user, an input device 40 that accepts various inputs by the user, and a display device 50 that displays an image generated by the processing device 30.
  • the TOF sensor unit 10 includes an illuminating unit 11 that irradiates light whose light intensity is periodically modulated to an object to be measured, and a light receiving unit 12 that receives backscattered light from the object to be measured. Have.
  • FIG. 2 is a diagram illustrating a state in which the TOF sensor unit, the angle distribution measurement illumination unit, and the angle distribution measurement light receiving unit illustrated in FIG. 1 are arranged with respect to an object to be measured.
  • the TOF sensor unit 10 emits light so as to be substantially orthogonal to the surface of the measurement object 2 and receives backscattered light in a direction substantially orthogonal to the surface of the measurement object 2.
  • the light irradiation angle and the light reception angle are deviated from the direction orthogonal to the surface of the object 2 to be measured. ing.
  • the illumination unit 11 includes a light source that is, for example, an LED (Light Emitting Diode), and irradiates light having a wavelength of 850 nm.
  • the illumination unit 11 emits light whose light intensity has been modulated so as to be a 30 MHz sine wave in accordance with control by the processing device 30.
  • the intensity modulation of the light emitted from the light source may be pulsed light or a rectangular wave.
  • the wavelength of the light irradiated from the illuminating unit 11 is not particularly limited, and may be variable or light having a wide wavelength band.
  • the light irradiated from the light source may be irradiated to the measurement object 2 via a diffusion lens, a condensing lens, an optical fiber, or the like.
  • the light receiving unit 12 includes a photodetector that is a CMOS (Complementary Metal Oxide Semiconductor) type TOF sensor chip.
  • the photodetector converts the received light into an electrical signal by photoelectric conversion.
  • the photodetector may be a CCD (Charge Coupled Devices) type TOF sensor chip.
  • the photodetector may receive the backscattered light from the measurement object 2 via a condenser lens, a wavelength filter that transmits only light of a predetermined wavelength, or the like.
  • the angle distribution measuring device 20 has an angle distribution measuring illumination unit 20A and an angle distribution measuring light receiving unit 20B.
  • the angle distribution measurement illumination unit 20 ⁇ / b> A has an illumination unit 21 and irradiates light to the object to be measured 2 with variable angles.
  • the angle distribution measurement light receiving unit 20 ⁇ / b> B has a light receiving unit 22 and receives backscattered light of the irradiated light in the measurement object 2 in a variable angle.
  • the illumination part 21 has a light source which is LED, for example, and irradiates light with a wavelength of 850 nm.
  • the wavelength of the light irradiated from the illumination part 21 is not specifically limited, It may be variable and may be light with a wide wavelength band.
  • the light irradiated from the light source may be irradiated to the measurement object 2 via a diffusion lens, a condensing lens, an optical fiber, or the like. Further, similarly to the illumination unit 11, the illumination unit 21 may irradiate light whose light intensity is periodically modulated.
  • the light receiving unit 22 may have the same configuration as the light receiving unit 12, but may be different in configuration and sensitivity from the light receiving unit 12.
  • the processing device 30 includes a control unit 31, a TOF measurement unit 32, an angle distribution measurement unit 33, an inverse analysis calculation unit 34, an image generation unit 35, and a storage unit 36.
  • the control unit 31 is realized by using a CPU (Central Processing Unit) or the like.
  • the control unit 31 controls the processing operation of each unit of the processing device 30.
  • the control unit 31 controls the operation of the processing device 30 by transferring instruction information and data to each component of the processing device 30.
  • the control unit 31 performs control so that the operations of the illumination unit 11 and the light receiving unit 12 and the operations of the illumination unit 21 and the light receiving unit 22 are synchronized.
  • the TOF measurement unit 32 uses the intensity of the light irradiated from the illumination unit 11 and the light received by the light receiving unit 12 as a first measurement value obtained by measuring isotropic backscattered light of the light irradiated to the measurement object 2.
  • the TOF information that is information on the phase shift of the modulation is measured. Specifically, the phase of the light emitted from the illumination unit 11 is acquired from the control unit 31, the phase of the light received from the light receiving unit 12 is acquired, and the acquired phase is compared to calculate the phase shift ⁇ .
  • the calculated TOF information is stored in the storage unit 36.
  • the first measurement value may be a value measured using at least one of a photoacoustic phenomenon and a crossed Nicols method.
  • FIG. 3 is a diagram illustrating a state in which the first measurement value is measured by the TOF sensor unit illustrated in FIG.
  • the TOF sensor unit 10 emits light so that the illumination unit 11 is orthogonal to the surface of the measurement object 2, and the light receiving unit 12 a and the light receiving unit 12 b emit backscattered light of the measurement object 2.
  • Light is received in a direction perpendicular to the surface of the measurement object 2.
  • the backscattered light from the shallow layer of the measurement object 2 detected by the light receiving unit 12a located in the vicinity of the illumination unit 11 is anisotropic backscattered light with many components in the direction orthogonal to the surface of the measurement object 2. It is.
  • the backscattered light from the deep layer of the measurement object 2 detected by the light receiving unit 12b disposed at a position farther from the illumination unit 11 than the light receiving unit 12a is isotropically back with components in each direction made uniform. Scattered light.
  • the TOF sensor unit 10 is suitable for detecting information from a deeper layer by detecting information on the object 2 to be measured based on the phase shift ⁇ .
  • the TOF sensor unit 10 may have a plurality of light receiving portions corresponding to the light receiving portions 12a and 12b, and one light receiving portion corresponds to the light receiving portions 12a and 12b. The structure which moves to a position and measures may be sufficient.
  • the angle distribution measurement unit 33 performs backscattering on the light irradiated from the illumination unit 21 as a second measurement value obtained by measuring backscattered light that is more anisotropic than the first measurement value of the light irradiated on the measurement object 2. Measure the angular distribution of light intensity. Specifically, under the control of the control unit 31, the light intensity is detected by changing the angle ⁇ 2 received by the light receiving unit 22 with respect to the light emitted from the illumination unit 21 at the angle ⁇ 1. Further, the angle ⁇ 1 irradiated from the illuminating unit 21 is changed, and the angle ⁇ 2 received by the light receiving unit 22 is similarly changed to repeat the measurement. The measured light intensity angle distribution is stored in the storage unit 36.
  • the second measured value may be a value measured using at least one of the OCT (Optical Coherence Tomography) method and the parallel Nicol method.
  • OCT Optical Coherence Tomography
  • parallel Nicol With orthogonal Nicol, it is possible to measure multiple scattered light within the object 2 to be measured.
  • parallel Nicol light that has been multiple scattered within the object 2 and reflected or scattered at the surface of the object 2 (single scattering). The measured light can be measured. Therefore, the second measured value may be a value measured by using both the orthogonal Nicol method and the parallel Nicol method.
  • FIG. 4 is a diagram illustrating a state in which the second measurement value is measured by the angle distribution measuring apparatus illustrated in FIG.
  • the illumination unit 21 irradiates the surface of the measurement object 2 with light at an angle ⁇ ⁇ b> 1, and the light receiving unit 22 a and the light receiving unit 22 b are backscattered of the measurement object 2.
  • Light is received at an angle ⁇ 2 with respect to the surface of the object 2 to be measured.
  • the backscattered light from the shallow layer of the DUT 2 detected by the light receiving unit 22a located in the vicinity of the illuminating unit 21 has a component in a specific direction including information on the light propagation path and the optical path length of the backscattered light.
  • the angle distribution measuring apparatus 20 can obtain shallow layer information with high accuracy, it is not suitable for obtaining deep layer information.
  • the inverse analysis calculation unit 34 reads the TOF information measured by the TOF measurement unit 32 from the storage unit 36 and estimates a first optical property value distribution that is an optical property value distribution inside the measurement object 2 by inverse analysis calculation. . Further, the inverse analysis calculation unit 34 reads the light intensity angle distribution of the backscattered light measured by the angle distribution measurement unit 33 from the storage unit 36, and performs inverse using at least a part of the first optical physical property value distribution as an initial value. A second optical property value distribution, which is a shallow optical property value distribution, is estimated from the first optical property value distribution inside the object to be measured 2 by analytical calculation.
  • the inverse analysis calculation unit 34 may cause a CPU or the like provided inside the optical property value distribution estimation apparatus 1 to perform calculations, but a CPU or the like provided outside the optical property value distribution estimation apparatus 1 An inverse analysis operation may be performed.
  • first optical property value distribution and the second optical property value distribution are, for example, scattering coefficient distributions. Further, the first optical property value distribution and the second optical property value distribution may be an absorption coefficient distribution, a scattering anisotropy distribution, or the like.
  • the first optical property value distribution which is the property value distribution
  • the second optical property value distribution which is the optical property value distribution of the second region (shallow layer) that is estimated as follows, is estimated.
  • the image generation unit 35 generates and outputs an image signal based on the optical characteristics obtained by the inverse analysis calculation unit 34 by calculation. Specifically, the image generation unit 35 generates an image signal in which the density, color, contrast, and the like are changed according to the value of the optical characteristic for each pixel obtained by the inverse analysis calculation unit 34 by calculation.
  • the image generation unit 35 may generate an image signal in which an image whose color or the like is changed according to the value of optical characteristics is superimposed on a 2D or 3D image captured by various cameras, an ultrasonic image, or the like. Further, the image generation unit 35 may generate an image in which the inverse analysis calculation unit 34 expresses the optical characteristics with numerical values, graphs, or the like by calculation.
  • the storage unit 36 stores various programs for operating the optical property value distribution estimation apparatus 1, data including various parameters necessary for the operation of the optical property value distribution estimation apparatus 1, and the like. Further, the storage unit 36 stores the first measurement value based on the measurement result of the TOF measurement unit 32 and also stores the second measurement value based on the measurement result of the angle distribution measurement unit 33.
  • the storage unit 36 stores various programs including an operation program for executing the operation method of the optical property value distribution estimation apparatus 1.
  • the operation program can be stored in a computer-readable storage medium such as a hard disk, a flash memory, a CD-ROM, a DVD-ROM, or a flexible disk and widely distributed.
  • the various programs described above can also be obtained by downloading via a communication network.
  • the communication network here is realized by, for example, an existing public line network, a LAN (Local Area Network), a WAN (Wide Area Network) or the like, and may be wired or wireless.
  • the storage unit 36 having the above configuration is realized using a ROM (Read Only Memory) in which various programs are installed in advance, and a RAM (Random Access Memory) that stores calculation parameters and data of each process. .
  • ROM Read Only Memory
  • RAM Random Access Memory
  • the input device 40 is realized by using an operation device such as a mouse, a keyboard, and a touch panel, and receives input of various instruction information of the optical property value distribution estimation device 1 and outputs it to the control unit 31.
  • an operation device such as a mouse, a keyboard, and a touch panel
  • the display device 50 is configured by using a display or the like using liquid crystal or organic EL (Electro Luminescence).
  • the display device 50 displays an image based on the image signal output from the image generation unit 35.
  • FIG. 5 is a flowchart showing an operation of measuring an object to be measured by the optical property value distribution estimating apparatus shown in FIG. As shown in FIG. 5, first, intensity-modulated light is irradiated onto the measurement object 2 from the illumination unit 11 of the TOF sensor unit 10 (step S1). Then, the light receiving unit 12 of the TOF sensor unit 10 receives the backscattered light from the measurement object 2 (step S2).
  • the TOF measurement unit 32 of the processing device 30 is information on the phase shift ⁇ of intensity modulation between the light emitted from the illumination unit 11 and the light received by the light reception unit 12 based on the measurement value of the light reception unit 12.
  • the TOF information is measured (step S3).
  • the inverse analysis calculation unit 34 estimates the first optical property value distribution inside the measurement object 2 by the reverse analysis calculation based on the TOF information (step S4: first estimation step).
  • the object to be measured 2 is irradiated with light at a predetermined irradiation angle from the illumination unit 21 of the angle distribution measurement illumination unit 20A (step S5).
  • the light receiving unit 22 of the angle distribution measurement light receiving unit 20B receives the backscattered light from the measurement object 2 at a predetermined light receiving angle (step S6).
  • control unit 31 determines whether or not the measurement has been completed at all light receiving angles (step S7).
  • step S7 determines that the measurement has not been completed at all the light reception angles (step S7: No)
  • step S8 the angle at which the light reception unit 22 receives light is changed under the control of the control unit 31 (step S8). The measurement in step S5 and step S6 is repeated.
  • step S7 determines whether the measurement is completed at all the light receiving angles.
  • step S9 determines whether the measurement is completed at all the irradiation angles. ).
  • step S9 determines that the measurement has not been completed at all the irradiation angles.
  • step S10 under the control of the control unit 31, the angle at which the illumination unit 21 emits light is changed (step). S10), the measurements in steps S5 to S7 are repeated.
  • step S9 determines that the measurement has been completed at all irradiation angles
  • step S9: Yes the angle distribution measurement unit 33 of the processing device 30 is based on the measurement value of the light receiving unit 22, and the illumination unit The angular distribution of the light intensity of the backscattered light with respect to the light irradiated from 21 is measured (step S11).
  • the inverse analysis calculation unit 34 performs an inverse analysis calculation using the first optical property value distribution as an initial value, thereby measuring the object 2 to be measured.
  • a second optical property value distribution which is a shallow optical property value distribution, is estimated from the internal first optical property value distribution (step S12: second estimation step).
  • step S1 to step S3 the measurement in step S1 to step S3 is performed, the calculation in step S4 is performed, the measurement in step S5 to step S11 is performed, and the calculation in step S12 is performed is described.
  • the measurement in step S1 to step S3 and the measurement in step S5 to step S11 may be performed in order or in parallel, and then the calculation in step S4 and the calculation in step S12 may be performed.
  • the inverse analysis calculation unit 34 estimates the first optical property value distribution inside the measurement object 2 by the inverse analysis calculation based on the TOF information. . Thereafter, the inverse analysis calculation unit 34 estimates the second optical property value distribution, which is the optical property value of the shallow layer, by inverse analysis calculation using the first optical property value distribution as an initial value based on the light intensity angle distribution. To do.
  • the light propagation component in the shallow layer consists of a direct propagation component due to the light from the illumination light and an indirect propagation component in which the return light due to the backscattering of the light reaching the deep layer reilluminates the shallow layer. Can think. For this reason, in order to estimate the optical property value of the shallow region with high accuracy, it is necessary to roughly know the intensity of the backscattered light from the deep region.
  • the propagating component is more strongly affected by the optical property value distribution in the deep layer region than in the shallow layer region.
  • FIG. 6 is a diagram showing light propagation in the deep layer and the shallow layer.
  • the deep layer region has a larger volume (area in the cross section of FIG. 6) than the shallow layer region. For this reason, the deep layer region has a greater influence on the light intensity of the backscattered light received by the light receiving unit 12 as a whole region than the shallow layer region. That is, in the process of performing the estimation calculation, the deep layer region greatly contributes to the error function.
  • the influence of the backscattered light received by the light receiving unit 12 on the light intensity by changing the physical property value of each mesh is Is smaller than the shallow layer.
  • the sensitivity of the change in the light intensity of the backscattered light to the change in the physical property value of each mesh is higher in the shallow layer than in the deep layer. That is, in the process of performing the estimation calculation, when viewed relatively, the physical property value parameter of each mesh in the shallow layer changes rapidly, and the change in the physical property value parameter of each mesh in the deep layer slows down. For this reason, if the physical property value parameters of the shallow layer and the deep layer are estimated and calculated at the same time (or the estimation calculation is performed first from the shallow layer), only the physical property parameter of the shallow layer changes, and the deep layer cannot be estimated. In addition, since the accuracy of calculating the amount of backscattered light from the deep layer to the shallow layer is deteriorated, the estimation accuracy of the shallow layer is also lowered.
  • the optical property value distribution estimation apparatus 1 first estimates the first optical property value in a wide range from the shallow layer to the deep layer based on the TOF information.
  • the reason why the TOF information is used for estimating the optical property value in the deep layer region and the light intensity angle distribution is used for estimating the optical property value in the shallow layer region will be described. Since the backscattered light returning from the shallow region to the surface layer has a small number of scattering, information on the light propagation path and the optical path length remains in the information on the light intensity angle distribution. On the other hand, the backscattered light that has reached the deep layer region and returned to the surface layer has a large number of scattering, the light intensity angular distribution becomes isotropic, and information on the light propagation path and the optical path length is lost. On the other hand, information on the light propagation path and the optical path length remains in the TOF information.
  • TOF information to estimate the optical property value distribution in the deep layer region and to use angle distribution information of the backscattered light to estimate the optical property value distribution in the shallow layer region.
  • Various known methods may be used for detecting the backscattered light from the shallow layer region and the deep layer region.
  • FIG. 7 is a diagram illustrating an object to be measured that is measured by the optical property value distribution estimating apparatus according to the embodiment.
  • the optical property value distribution estimation apparatus 1 is used to measure the optical property value distribution of the measurement object 2 having different optical property values between the first layer and the second layer.
  • the thickness of the first layer is 2 mm
  • the thickness of the entire object to be measured 2 is 7 mm
  • the width is 40 mm.
  • FIG. 8 is a diagram showing optical property values of the object to be measured shown in FIG.
  • the absorption coefficient ⁇ a [1 / mm] and the scattering anisotropy g are equal, but the scattering coefficient ⁇ s [1 / mm] is different.
  • FIG. 9 is a diagram showing a measurement result by the TOF measurement unit.
  • the distance l1 [mm] (see FIG. 2) between the illumination unit 11 and the light receiving unit 12 is increased, the intensity of the light emitted from the illumination unit 11 and the light received by the light receiving unit 12 is increased.
  • the amount [ps] of the modulation phase shift ⁇ increases.
  • the phase shift ⁇ may be measured by other known techniques.
  • FIG. 10 is a diagram illustrating a state in which the object to be measured illustrated in FIG. 7 is divided into meshes. As shown in FIG. 10, the measurement object 2 is divided into meshes.
  • FIG. 11 is an enlarged view of the mesh shown in FIG. As shown in FIG. 11, each mesh is further divided into 16 regions divided by an angle ⁇ (22.5 °).
  • FIG. 12 is a diagram illustrating a state in which an illumination unit and a light receiving unit are arranged on the measurement object illustrated in FIG.
  • the illumination unit 11 disposed in the center of the measurement object 2 irradiates light at an irradiation angle of 90 ° so as to be orthogonal to the surface of the measurement object 2.
  • the light irradiated from the illumination part 11 shall spread at +/- 11.25 degrees.
  • the light receiving unit 12 disposed on the entire surface of the measurement object 2 receives light from all directions at a light reception angle of 90 ° to ⁇ 90 ° orthogonal to the surface of the measurement object 2.
  • the optical property value distribution may be estimated using a diffusion equation instead of the light transport equation.
  • numerical calculation is performed in advance using a light propagation model assuming an optical physical property value distribution, and the optical physical property value distribution and a calculated value of backscattered light corresponding thereto are set as a lookup table, and simulated annealing is performed.
  • the optical property value distribution may be estimated using an algorithm such as a method.
  • FIG. 13 is a diagram showing an estimation result by the first estimation step.
  • the estimation result of the scattering coefficient is calculated by performing an inverse analysis calculation assuming that the initial value of the scattering coefficient is 1. This estimation result is substantially equal to the actual value in the deep layer (second layer).
  • FIG. 14 is a diagram for explaining an initial value setting method in the second estimation step.
  • the optical property values of the regions other than the region A are replaced by the estimation result of the columnar region A immediately below the illumination unit 11 to obtain the initial value of the second step.
  • the second optical property value distribution is the optical property of the vertical region A in the direction perpendicular to the surface of the measurement object 2 along the light irradiated to the measurement object 2 in the first optical property value distribution.
  • a value obtained by replacing the optical property values of the entire region of the object to be measured 2 is estimated as an initial value. This is because the region A is directly under the illumination unit 11 and the optical property value is estimated with the highest accuracy.
  • an inverse analysis calculation using a part of the first optical property value distribution as an initial value may be performed.
  • the optical property value used in the second estimation step is preferentially calculated, the calculation of other optical property values in the first estimation step, the second estimation step, May be performed in parallel.
  • an optical property value of a region other than this region may be replaced with the initial value of the second step based on the estimation result of the columnar region along the light irradiated from the illumination unit 11.
  • the illumination unit 11 and the light receiving unit 12 include not only one but also a plurality of light sources and photodetectors.
  • the optical property values of the regions other than these regions are replaced, and the initial stage of the second estimation step It can be a value.
  • any one of the estimation results of a plurality of columnar regions may be used, or a value calculated by averaging or weighting the estimation results of each columnar region
  • values calculated by various known methods may be used.
  • 15 to 18 are diagrams showing measurement results obtained by the angle distribution measurement unit. 15 to 18, the distance l2 (see FIG. 2) between the illumination unit 21 and the light receiving unit 22 is changed from 0 mm to 1.0 mm, and the light receiving angle of the light receiving unit 22 is 0 ° to 60 ° at each distance l2. Represents a part of the measured value when changed to. As described above, when the distance l2 between the illumination unit 21 and the light receiving unit 22 and the light receiving angle are changed, the angular distribution of the light intensity can be measured.
  • the light receiving angle is the central angle of light received by the light receiving unit 22, and the light receiving unit 22 receives light having a light receiving angle of ⁇ 11.25 °.
  • the calculation model shown in FIG. 12 is used to calculate the state of light propagating through the measurement object 2, and the optical property value is calculated.
  • the light receiving unit 22 arranged on the entire surface of the measurement object 2 performs calculation by changing the light receiving angle from 0 ° to 60 ° so as to correspond to the measured value.
  • FIG. 19 is a diagram illustrating an estimation result obtained in the second estimation step. As shown in FIG. 19, the estimation result of the scattering coefficient substantially equal to the actual value could be obtained by estimating the optical property value using the value of the first estimation step as the initial value.
  • FIG. 20 is a diagram illustrating an example of a method for estimating the thickness of the first layer. As shown in FIG. 20, when the approximate equivalent scattering coefficient ⁇ s ′ round of the first layer of the measurement object 2 is known, the thickness of the first layer can be calculated by, for example, 1 / ⁇ s ′ round. it can.
  • the depth of the region having a small optical property value in the surface layer region of the object to be measured 2 is predicted, and the angle at which the light to be measured 2 is irradiated according to the depth, and the rear from the object 2 to be measured
  • the optical property value can be estimated with high accuracy in a shorter time. That is, it is preferable that the second measurement value is measured by a measurement method determined in accordance with an optical property value (scattering coefficient or the like) representing the known scattering intensity.
  • the configuration in which the TOF measurement unit 32 calculates the intensity modulation phase shift between the light emitted from the illumination unit 11 and the light received by the light receiving unit 12 as TOF information has been described. Not limited to.
  • the TOF measurement unit 32 is configured as TOF information. The time delay of the light received by the light receiving unit 12 with respect to the pulsed light emitted from the illumination unit may be calculated.
  • SYMBOLS 1 Optical physical property value distribution estimation apparatus 2 To-be-measured object 10 TOF sensor unit 11, 21 Illumination part 12, 12a, 12b, 22, 22a, 22b Light-receiving part 20 Angle distribution measurement apparatus 20A Angle distribution measurement illumination unit 20B Angle distribution measurement light-receiving unit DESCRIPTION OF SYMBOLS 30 Processing apparatus 31 Control part 32 TOF measurement part 33 Angle distribution measurement part 34 Inverse analysis calculation part 35 Image generation part 36 Storage part 40 Input device 50 Display apparatus

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Pathology (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Engineering & Computer Science (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Optics & Photonics (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置が行う光学物性値分布の推定方法であって、前記被計測物に照射した光の等方的な後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光学物性値分布である第１の光学物性値分布を推定する第１の推定工程と、前記被計測物に照射した光の前記第１の計測値より非等方的な後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の光学物性値分布である第２の光学物性値分布を推定する第２の推定工程と、を含む。これにより、生体組織の浅層から深層に至る広い範囲に渡って光学物性値分布を高精度に推定することができる光学物性値分布の推定方法を提供する。

Description

光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置
　本発明は、光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置に関する。
　従来、生体組織内部における光学物性値の分布を非侵襲に推定する方法として、反射型拡散光トモグラフィー法を用いた方法が知られている（例えば、特許文献１～３参照）。この方法では、はじめに、被計測物に光を照射して、被計測物からの後方散乱光を受光して、後方散乱光の計測値を取得する。
　続いて、被計測物内部の光学物性値分布を仮定して光伝播モデルを用いて数値計算を行い、この被計測物による後方散乱光を予測値として算出する。具体的には、被計測物モデルをメッシュ（又はボクセル）に区切り、それぞれのメッシュに対して散乱係数と吸収係数とをパラメータとして設定し、被計測物内部を光が伝播する様子を光輸送方程式や光拡散方程式により算出する。そして、測定値と予測値とを比較し、一致度が所定値以上であるか否かを判定する。一致度が所定値より低い場合には、被計測物内部の光学物性値分布を再度仮定し直して予測値を算出し、一致度が所定値以上になるまでこの計算を繰り返す。この繰り返し計算は逆解析演算とも呼ばれる。なお、計測値と予測値との差から、誤差関数を設定することにより一致度を評価することができる。この誤差関数が小さくなるように、被計測物内部の光学物性値分布を再度仮定し直すことにより、測定値と予測値との差を小さくすることができる。上述した逆解析演算には、公知の様々な最適化アルゴリズムを適用することができる。
　特許文献１には、被計測物に対して斜めに照明を照射し、照明方向に対向する斜めの方向において被計測物からの後方散乱光を受光することにより、被計測物の浅層の光学物性値分布を高精度に推定できることが開示されている。
　特許文献２には、複数の照明方向及び複数の受光方向において計測を行うことにより、被計測物の浅層の光学物性値分布をさらに高精度に推定できることが開示されている。
　特許文献３には、生体の深層に到達した光の後方散乱光は、散乱回数が多いために等方的になっていることが開示されている。換言すると、生体の深層に到達した光の後方散乱光は、光強度の角度依存性が低い。一方で、被計測物からの後方散乱光のＴＯＦ（Ｔｉｍｅ　Ｏｆ　Ｆｌｉｇｈｔ）情報である位相ずれは、深層に到達した光ほど大きくなる傾向があることにより、ＴＯＦ情報を用いて被計測物の深層の光学物性値分布を高精度に推定できることが開示されている。
特表２００５－５１３４９１号公報 国際公開第２００２／０６９７９２号 特開２００６－２００９４３号公報
　ところで、医療用イメージング機器や医療用診断装置の開発においては、生体組織の浅層から深層に至る広い範囲に渡って光学物性値分布（散乱係数、吸収係数等の分布）を知りたいという要望がある。
　特許文献１～３には、生体組織の浅層、深層のそれぞれの推定に適した計測方法が開示されているものの、浅層から深層に至る広い範囲に渡って光学物性値分布を高精度に推定することができないという課題があった。
　本発明は、上記に鑑みてなされたものであって、生体組織の浅層から深層に至る広い範囲に渡って光学物性値分布を高精度に推定することができる光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置を提供することを目的とする。
　上述した課題を解決し、目的を達成するために、本発明の一態様に係る光学物性値分布の推定方法は、被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置が行う光学物性値分布の推定方法であって、前記被計測物に照射した光の等方的な後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光学物性値分布である第１の光学物性値分布を推定する第１の推定工程と、前記被計測物に照射した光の前記第１の計測値より非等方的な後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の光学物性値分布である第２の光学物性値分布を推定する第２の推定工程と、を含むことを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、第２の光学物性値分布は、前記被計測物内部の前記第１の光学物性値より浅層の光学物性値分布であることを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、前記第２の推定工程は、前記第１の光学物性値分布のうち、前記被計測物に照射した光に沿って前記被計測物の表面に垂直な方向の垂直領域の光学物性値により、前記被計測物の全領域の光学物性値を置換した値を初期値として前記第２の光学物性値分布を推定することを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、前記第１の計測値は、前記被計測物にパルス光又は周期的に強度を変調させた光を照射して計測したＴＯＦ情報を含む値であり、前記第１の推定工程は、前記ＴＯＦ情報に基づいて前記第１の光学物性値分布を推定することを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、前記第２の計測値は、前記被計測物に光を照射する角度、又は前記被計測物からの後方散乱光を受光する角度の少なくとも一方を、前記被計測物の表面に対して鋭角に角度をつけて計測した値であり、前記第２の推定工程は、前記第２の計測値の光強度角度分布に基づいて前記第２の光学物性値分布を推定することを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、前記第１の計測値は、光音響現象、直交ニコル法の少なくともいずれか１つを用いて計測した値であることを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、前記第２の計測値は、ＯＣＴ法、平行ニコル法の少なくともいずれか１つを用いて計測した値であることを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、前記第２の計測値は、既知である散乱の強さを表す光学物性値に応じて定められた計測方法により計測した値であることを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、前記計測方法は、前記被計測物の表層領域における前記光学物性値が小さい領域の深さを予測し、該深さに応じて前記被計測物に光を照射する角度、及び前記被計測物からの後方散乱光を受光する角度を定め、該定められた角度で前記第２の計測値を計測することを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、前記深さは、前記被計測物の概略の等価散乱係数μｓ’ｒｏｕｇｈを用いて、１／μｓ’ｒｏｕｇｈにより算出されることを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置が行う光学物性値分布の推定方法であって、前記被計測物に照射した光の後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光エネルギーの損失が大きいと推測される領域を含む第１の領域の光学物性値分布である第１の光学物性値分布を推定する第１の推定工程と、前記被計測物に照射した光の後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の前記第１の領域より光エネルギーの損失が小さいと推測される第２の領域の光学物性値分布である第２の光学物性値分布を推定する第２の推定工程と、を含むことを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定方法は、被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置が行う光学物性値分布の推定方法であって、前記被計測物の第１の領域に照射した光の後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光学物性値分布である第１の光学物性値分布を推定する第１の推定工程と、前記第１の領域に包含される第２の領域に照射した光の後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の前記第１の光学物性値分布より浅層の光学物性値分布である第２の光学物性値分布を推定する第２の推定工程と、を含むことを特徴とする。
　また、本発明の一態様に係る光学物性値分布の推定プログラムは、被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置が行う光学物性値分布の推定プログラムであって、前記被計測物に照射した光の等方的な後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光学物性値分布である第１の光学物性値分布を推定する第１の推定ステップと、前記被計測物に照射した光の前記第１の計測値より非等方的な後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の前記第１の光学物性値分布より浅層の光学物性値分布である第２の光学物性値分布を推定する第２の推定ステップと、をコンピュータに実行させることを特徴とする。
　また、本発明の一態様に係る光学物性値分布推定装置は、被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置であって、前記被計測物に照射した光の等方的な後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光学物性値分布である第１の光学物性値分布を推定し、前記被計測物に照射した光の前記第１の計測値より非等方的な後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の前記第１の光学物性値分布より浅層の光学物性値分布である第２の光学物性値分布を推定する逆解析演算部を備えることを特徴とする。
　本発明によれば、生体組織の浅層から深層に至る広い範囲に渡って光学物性値分布を高精度に推定することができる光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置を実現することができる。
図１は、実施の形態１に係る光学物性値分布推定装置の構成を示すブロック図である。 図２は、図１に示すＴＯＦセンサユニット、角度分布計測照明ユニット、及び角度分布計測受光ユニットを被計測物に対して配置した様子を表す図である。 図３は、図２に示すＴＯＦセンサユニットにより第１の計測値を計測する様子を表す図である。 図４は、図２に示す角度分布計測装置により第２の計測値を計測する様子を表す図である。 図５は、図１に示す光学物性値分布推定装置により被計測物を計測する動作を示すフローチャートである。 図６は、深層及び浅層における光の伝播を表す図である。 図７は、実施例に係る光学物性値分布推定装置により計測を行う被計測物を表す図である。 図８は、図７に示す被計測物の光学物性値を表す図である。 図９は、ＴＯＦ計測部による計測結果を表す図である。 図１０は、図７に示す被計測物をメッシュ状に分割する様子を表す図である。 図１１は、図１０に示すメッシュの拡大図である。 図１２は、図１０に示す被計測物に照明部と受光部とを配置した様子を表す図である。 図１３は、第１の推定工程による推定結果を表す図である。 図１４は、第２の推定工程による初期値の設定方法を説明するための図である。 図１５は、角度分布計測部による計測結果を表す図である。 図１６は、角度分布計測部による計測結果を表す図である。 図１７は、角度分布計測部による計測結果を表す図である。 図１８は、角度分布計測部による計測結果を表す図である。 図１９は、第２の推定工程による推定結果を表す図である。 図２０は、第１層の厚さを推定する方法の一例を表す図である。
　以下に、図面を参照して本発明に係る光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置の実施の形態を説明する。なお、これらの実施の形態により本発明が限定されるものではない。以下の実施の形態においては、計測値としてＴＯＦ情報及び光強度角度分布を用いて光学物性値分布を推定する光学物性値分布の推定方法を例示して説明するが、本発明は、他の計測値を用いる光学物性値分布の推定方法にも適用することができる。
　また、図面の記載において、同一又は対応する要素には適宜同一の符号を付している。また、図面は模式的なものであり、各要素の寸法の関係、各要素の比率等は、現実と異なる場合があることに留意する必要がある。図面の相互間においても、互いの寸法の関係や比率が異なる部分が含まれている場合がある。
（実施の形態１）
　図１は、実施の形態１に係る光学物性値分布推定装置の構成を示すブロック図である。本実施の形態１に係る光学物性値分布推定装置１は、被計測物の内部の光学物性値分布を推定する装置である。光学物性値分布推定装置１は、図１に示すように、被計測物に対して周期的に強度を変調させた光を照射するとともに、被計測物における後方散乱光を受光するＴＯＦセンサユニット１０と、被計測物に対して角度を可変に光を照射するとともに、被計測物における後方散乱光を、角度を可変に受光する角度分布計測装置２０と、ＴＯＦセンサユニット１０及び角度分布計測装置２０が取得した信号を処理する処理装置３０と、ユーザによる各種入力を受け付ける入力装置４０と、処理装置３０が生成した画像を表示する表示装置５０と、を備える。
　ＴＯＦセンサユニット１０は、被計測物に対して周期的に光強度を変調させた光を照射する照明部１１と、照射した光の被計測物における後方散乱光を受光する受光部１２と、を有する。
　図２は、図１に示すＴＯＦセンサユニット、角度分布計測照明ユニット、及び角度分布計測受光ユニットを被計測物に対して配置した様子を表す図である。ＴＯＦセンサユニット１０は、被計測物２の表面に対して略直交するように光を照射し、被計測物２の表面に対して略直交する方向の後方散乱光を受光する。なお、図２では、記載の都合上、光を照射する角度及び受光する角度が被計測物２の表面に対して直交する方向からずれているが、実際には、略直交するように配置されている。
　照明部１１は、例えば、ＬＥＤ（Ｌｉｇｈｔ　Ｅｍｉｔｔｉｎｇ　Ｄｉｏｄｅ）である光源を有し、波長８５０ｎｍの光を照射する。照明部１１は、処理装置３０による制御に応じて、３０ＭＨｚの正弦波となるように光強度が変調された光を照射する。ただし、光源から照射する光の強度変調は、パルス光や矩形波としてもよい。また、照明部１１から照射する光の波長は、特に限定されず、可変であってもよいし、波長帯域の広い光であってもよい。また、光源から照射された光は、拡散レンズや集光レンズ、光ファイバ等を介して被計測物２に照射されてもよい。
　受光部１２は、ＣＭＯＳ（Ｃｏｍｐｌｅｍｅｎｔａｒｙ　Ｍｅｔａｌ　Ｏｘｉｄｅ　Ｓｅｍｉｃｏｎｄｕｃｔｏｒ）型のＴＯＦセンサチップである光検出器を有する。光検出器は、受光した光を光電変換により電気信号に変換する。ただし、光検出器は、ＣＣＤ（Ｃｈａｒｇｅ　Ｃｏｕｐｌｅｄ　Ｄｅｖｉｃｅｓ）型のＴＯＦセンサチップであってもよい。また、光検出器は、被計測物２における後方散乱光を集光レンズや、所定の波長の光のみを透過させる波長フィルタ等を介して受光してもよい。
　角度分布計測装置２０は、角度分布計測照明ユニット２０Ａと、角度分布計測受光ユニット２０Ｂと、を有する。角度分布計測照明ユニット２０Ａは、照明部２１を有し、被計測物２に対して角度を可変に光を照射する。角度分布計測受光ユニット２０Ｂは、受光部２２を有し、照射した光の被計測物２における後方散乱光を角度を可変に受光する。
　照明部２１は、例えば、ＬＥＤである光源を有し、波長８５０ｎｍの光を照射する。また、照明部２１から照射する光の波長は、特に限定されず、可変であってもよいし、波長帯域の広い光であってもよい。また、光源から照射された光は、拡散レンズや集光レンズ、光ファイバ等を介して被計測物２に照射されてもよい。また、照明部２１は、照明部１１と同様に、周期的に光強度が変調された光を照射してもよい。
　受光部２２は、受光部１２と同様の構成であってもよいが、受光部１２と構成や感度が異なっていてもよい。
　処理装置３０は、制御部３１と、ＴＯＦ計測部３２と、角度分布計測部３３と、逆解析演算部３４と、画像生成部３５と、記憶部３６と、を有する。
　制御部３１は、ＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）等を用いて実現される。制御部３１は、処理装置３０の各部の処理動作を制御する。制御部３１は、処理装置３０の各構成に対する指示情報やデータの転送等を行うことによって、処理装置３０の動作を制御する。また、制御部３１は、照明部１１と受光部１２との動作、及び照明部２１と受光部２２との動作が、それぞれ同期するように制御を行う。
　ＴＯＦ計測部３２は、被計測物２に照射した光の等方的な後方散乱光を計測した第１の計測値として、照明部１１から照射する光と受光部１２が受光した光との強度変調の位相ずれの情報であるＴＯＦ情報を計測する。具体的には、制御部３１から照明部１１が照射する光の位相を取得し、受光部１２から受光した光の位相を取得し、取得した位相を比較して位相ずれφを算出する。算出されたＴＯＦ情報は、記憶部３６に記憶される。なお、第１の計測値は、光音響現象、直交ニコル法の少なくともいずれか１つを用いて計測した値であってもよい。
　図３は、図２に示すＴＯＦセンサユニットにより第１の計測値を計測する様子を表す図である。図３に示すように、ＴＯＦセンサユニット１０は、照明部１１が被計測物２の表面に直交するように光を照射し、受光部１２ａ及び受光部１２ｂが被計測物２の後方散乱光を被計測物２の表面に直交する方向において受光する。照明部１１の近傍に位置する受光部１２ａが検出する被計測物２の浅層からの後方散乱光は、被計測物２の表面に直交する方向の成分が多い非等方的な後方散乱光である。一方、受光部１２ａより照明部１１から離れた位置に配置された受光部１２ｂが検出する被計測物２の深層からの後方散乱光は、各方向の成分が均一にされた等方的な後方散乱光である。ＴＯＦセンサユニット１０は、位相ずれφにより被計測物２の情報を検出することにより、より深層からの情報を検出するのに適している。なお、ＴＯＦセンサユニット１０は、受光部１２ａ及び受光部１２ｂに対応する複数の受光部を有している構成であってもよいし、１つの受光部を受光部１２ａ及び受光部１２ｂに対応する位置に移動させて計測する構成であってもよい。
　角度分布計測部３３は、被計測物２に照射した光の第１の計測値より非等方的な後方散乱光を計測した第２の計測値として、照明部２１から照射する光に対する後方散乱光の光強度の角度分布を計測する。具体的には、制御部３１の制御のもと、照明部２１から角度θ１で照射する光に対して、受光部２２が受光する角度θ２を変えて光強度を検出する。さらに、照明部２１から照射する角度θ１を変え、同様に受光部２２が受光する角度θ２を変えて計測を繰り返す。計測された光強度角度分布は、記憶部３６に記憶される。なお、角度θ１又は角度θ２のいずれか一方のみを変化させてもよい。また、角度θ１又は角度θ２のいずれか一方が鋭角な角度（０°より大きく９０°より小さい）であり、他方が９０°であってもよい。なお、第２の計測値は、ＯＣＴ（Ｏｐｔｉｃａｌ　Ｃｏｈｅｒｅｎｃｅ　Ｔｏｍｏｇｒａｐｈｙ）法、平行ニコル法の少なくともいずれか１つを用いて計測した値であってもよい。直交ニコルでは、被計測物２内で多重散乱した光を計測することができ、平行ニコルでは、被計測物２内で多重散乱した光と被計測物２の表面で反射又は散乱（単散乱）された光とを計測することができる。従って、第２の計測値は、直交ニコル法と平行ニコル法とを併用して計測した値であってもよい。
　図４は、図２に示す角度分布計測装置により第２の計測値を計測する様子を表す図である。図４に示すように、角度分布計測装置２０は、照明部２１が被計測物２の表面に対して角度θ１で光を照射し、受光部２２ａ及び受光部２２ｂが被計測物２の後方散乱光を被計測物２の表面に対して角度θ２で受光する。照明部２１の近傍に位置する受光部２２ａが検出する被計測物２の浅層からの後方散乱光は、後方散乱光の光伝播経路や光路長の情報を含んでいる特定の方向の成分が多い非等方的な後方散乱光であり、後方散乱光の光強度が角度によって異なる。一方、受光部２２ａより照明部２１から離れた位置に配置された受光部２２ｂが検出する被計測物２の深層からの後方散乱光は、各方向の成分が均一にされた等方的な後方散乱光であり、後方散乱光の光強度がどの角度でも等しくなる。従って、角度分布計測装置２０は、浅層の情報を高精度に得ることができるものの、深層の情報を得るのには適していない。
　逆解析演算部３４は、ＴＯＦ計測部３２が計測したＴＯＦ情報を記憶部３６から読み出して、逆解析演算により被計測物２内部の光学物性値分布である第１の光学物性値分布を推定する。さらに、逆解析演算部３４は、角度分布計測部３３が計測した後方散乱光の光強度角度分布を記憶部３６から読み出して、第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、被計測物２内部の第１の光学物性値分布より浅層の光学物性値分布である第２の光学物性値分布を推定する。なお、逆解析演算部３４は、光学物性値分布推定装置１の内部に設けられたＣＰＵ等に演算を行わせてもよいが、光学物性値分布推定装置１の外部に設けられたＣＰＵ等に逆解析演算を行わせてもよい。
　ここで、第１の光学物性値分布及び第２の光学物性値分布は、例えば散乱係数の分布である。また、第１の光学物性値分布及び第２の光学物性値分布は、吸収係数の分布、又は散乱異方性の分布等であってもよい。
　生体組織（散乱吸収体）の内部を光が伝播する過程で光は拡散していくため、散乱や吸収により光エネルギーが失われていく。光伝播が非等方的（非拡散的）であるときは、生体組織内における実効的な光伝播距離が相対的に短いため、光エネルギー損失も相対的に小さい。一方で、光伝播が等方的（拡散的）になったときは、生体組織内における実効的な光伝播距離が相対的に長くなるため、光エネルギー損失も相対的に大きくなる。そのため、ＴＯＦ計測部３２で取得する第１の計測値を用いて被計測物２内部の光エネルギー損失が大きいと推測される領域（深層）を含む第１の領域（浅層及び深層）の光学物性値分布である第１の光学物性値分布を推定するとともに、角度分布計測部３３で取得する第２の計測値を用いて被計測物２内部の第１の領域より光エネルギーの損失が少ないと推測される第２の領域（浅層）の光学物性値分布である第２の光学物性値分布を推定する。
　画像生成部３５は、逆解析演算部３４が演算により求めた光学特性に基づいて画像信号を生成し、出力する。具体的には、画像生成部３５は、逆解析演算部３４が演算により求めた画素ごとの光学特性の値に応じて、濃淡や色、コントラスト等を変化させた画像信号を生成する。画像生成部３５は、各種カメラで撮像した２Ｄや３Ｄの画像、超音波画像等に、光学特性の値に応じて色等を変化させた画像を重畳した画像信号を生成してもよい。また、画像生成部３５は、逆解析演算部３４が演算により光学特性を数値やグラフ等で表した画像を生成してもよい。
　記憶部３６は、光学物性値分布推定装置１を動作させるための各種プログラム、及び光学物性値分布推定装置１の動作に必要な各種パラメータ等を含むデータ等を記憶する。また、記憶部３６は、ＴＯＦ計測部３２の計測結果に基づいて、第１の計測値を記憶するととともに、角度分布計測部３３の計測結果に基づいて、第２の計測値を記憶する。
　また、記憶部３６は、光学物性値分布推定装置１の作動方法を実行するための作動プログラムを含む各種プログラムを記憶する。作動プログラムは、ハードディスク、フラッシュメモリ、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ、フレキシブルディスク等のコンピュータ読み取り可能な記憶媒体に記憶して広く流通させることも可能である。なお、上述した各種プログラムは、通信ネットワークを介してダウンロードすることによって取得することも可能である。ここでいう通信ネットワークは、例えば既存の公衆回線網、ＬＡＮ（Ｌｏｃａｌ　Ａｒｅａ　Ｎｅｔｗｏｒｋ）、ＷＡＮ（Ｗｉｄｅ　Ａｒｅａ　Ｎｅｔｗｏｒｋ）等によって実現されるものであり、有線、無線を問わない。
　以上の構成を有する記憶部３６は、各種プログラム等が予めインストールされたＲＯＭ（Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）、及び各処理の演算パラメータやデータ等を記憶するＲＡＭ（Ｒａｎｄｏｍ　Ａｃｃｅｓｓ　Ｍｅｍｏｒｙ）等を用いて実現される。
　入力装置４０は、マウス、キーボード及びタッチパネル等の操作デバイスを用いて実現され、光学物性値分布推定装置１の各種指示情報の入力を受け付け、制御部３１に出力する。
　表示装置５０は、液晶又は有機ＥＬ（Ｅｌｅｃｔｒｏ　Ｌｕｍｉｎｅｓｃｅｎｃｅ）を用いた表示ディスプレイ等を用いて構成される。表示装置５０は、画像生成部３５が出力した画像信号に基づいた画像を表示する。
　次に、光学物性値分布推定装置１の動作を説明する。図５は、図１に示す光学物性値分布推定装置により被計測物を計測する動作を示すフローチャートである。図５に示すように、まず、ＴＯＦセンサユニット１０の照明部１１から、被計測物２に強度変調光を照射する（ステップＳ１）。そして、ＴＯＦセンサユニット１０の受光部１２が、被計測物２からの後方散乱光を受光する（ステップＳ２）。
　続いて、処理装置３０のＴＯＦ計測部３２が、受光部１２の計測値に基づいて、照明部１１から照射する光と受光部１２が受光した光との強度変調の位相ずれφの情報であるＴＯＦ情報を計測する（ステップＳ３）。
　さらに、逆解析演算部３４が、ＴＯＦ情報に基づいて、逆解析演算により被計測物２内部の第１の光学物性値分布を推定する（ステップＳ４：第１の推定工程）。
　続いて、角度分布計測照明ユニット２０Ａの照明部２１から、被計測物２に所定の照射角度で光を照射する（ステップＳ５）。そして、角度分布計測受光ユニット２０Ｂの受光部２２が、被計測物２からの後方散乱光を所定の受光角度で受光する（ステップＳ６）。
　その後、制御部３１は、全ての受光角度において計測が終了したか否かを判定する（ステップＳ７）。制御部３１が、全ての受光角度において計測が終了していないと判定した場合（ステップＳ７：Ｎｏ）、制御部３１の制御のもと、受光部２２が受光する角度を変えて（ステップＳ８）、ステップＳ５、ステップＳ６の計測が繰り返し行われる。
　一方、制御部３１が、全ての受光角度において計測が終了したと判定した場合（ステップＳ７：Ｙｅｓ）、制御部３１は、全ての照射角度において計測が終了したか否かを判定する（ステップＳ９）。制御部３１が、全ての照射角度において計測が終了していないと判定した場合（ステップＳ９：Ｎｏ）、制御部３１の制御のもと、照明部２１が光を照射する角度を変えて（ステップＳ１０）、ステップＳ５～ステップＳ７の計測が繰り返し行われる。
　一方、制御部３１が、全ての照射角度において計測が終了したと判定した場合（ステップＳ９：Ｙｅｓ）、処理装置３０の角度分布計測部３３が、受光部２２の計測値に基づいて、照明部２１から照射する光に対する後方散乱光の光強度の角度分布を計測する（ステップＳ１１）。
　さらに、逆解析演算部３４が、角度分布計測部３３が計測した後方散乱光の光強度角度分布に基づいて、第１の光学物性値分布を初期値とした逆解析演算により、被計測物２内部の第１の光学物性値分布より浅層の光学物性値分布である第２の光学物性値分布を推定する（ステップＳ１２：第２の推定工程）。
　なお、上述したフローチャートでは、ステップＳ１～ステップＳ３の計測を行い、ステップＳ４の演算を行い、ステップＳ５～ステップＳ１１の計測を行い、ステップＳ１２の演算を行う例を説明したがこれに限られない。例えば、ステップＳ１～ステップＳ３の計測とステップＳ５～ステップＳ１１の計測とを順番に又は並行して行い、その後、ステップＳ４の演算とステップＳ１２の演算とを行ってもよい。
　以上説明したように、光学物性値分布推定装置１により、はじめに、逆解析演算部３４が、ＴＯＦ情報に基づいて、逆解析演算により被計測物２内部の第１の光学物性値分布を推定する。その後、逆解析演算部３４が、光強度角度分布に基づいて、第１の光学物性値分布を初期値とした逆解析演算により浅層の光学物性値である第２の光学物性値分布を推定する。
　ここで、はじめに深層からの等方的な後方散乱光（被計測物２に直交する方向の後方散乱光）を用いて光学物性値を推定する理由を説明する。浅層において光が伝播する成分は、照明光からの光による直接的な伝播成分と、深層まで到達した光の後方散乱による戻り光が浅層を再照明する間接的な伝播成分と、からなると考えることができる。このため、浅層領域の光学物性値を高精度に推定するためには、深層領域からの後方散乱光の強さが概略的にわかっている必要がある。
　また、浅層領域に比べて深層領域の方が相対的に散乱又は吸収が強い場合、又は浅層領域に比べて深層領域の方が厚い場合（消化器組織等）には、深層において光が伝播する成分は、浅層領域よりも深層領域の光学物性値分布の影響を強く受けると考えられる。
　図６は、深層及び浅層における光の伝播を表す図である。図６に示すように、生体において、一般に深層領域は、浅層領域よりも体積（図６の断面における面積）が大きい。そのため、浅層領域より深層領域の方が、領域全体として受光部１２が受光する後方散乱光の光強度に対して与える影響が大きい。すなわち、推定計算を行う過程において、深層領域の方が誤差関数に大きく寄与する。一方で、深層まで光が到達する過程で光エネルギーは減少しているので、各メッシュの物性値が変化することにより、受光部１２が受光する後方散乱光の光強度に与える影響は、深層の方が浅層より小さい。換言すると、浅層では深層よりも、各メッシュの物性値の変化に対する後方散乱光の光強度の変化の感度が高い。すなわち、推定計算を行う過程では、相対的にみると、浅層の各メッシュの物性値パラメータは変化が敏速であり、深層の各メッシュの物性値パラメータは変化が緩慢になる。このため、浅層と深層との物性値パラメータを同時に推定計算する（又は浅層から先に推定計算する）と、浅層の物性値パラメータしか変化しないこととなり、深層の推定ができなくなるばかりでなく、深層から浅層方向への後方散乱光の量を計算する精度が悪くなることから浅層の推定精度も低下してしまう。
　以上の理由により、物性値パラメータの変化が緩慢である深層の光学物性値を先に推定することが好ましい。従って、光学物性値分布推定装置１では、はじめにＴＯＦ情報により浅層から深層に至る広い範囲の第１の光学物性値を推定する。
　次に、深層領域の光学物性値の推定にＴＯＦ情報、浅層領域の光学物性値の推定に光強度角度分布をそれぞれ用いる理由について説明する。浅層領域から表層へ戻ってきた後方散乱光は散乱回数が少ないために、光強度角度分布の情報には光伝播経路や光路長の情報が残っている。一方で、深層領域まで到達して表層へ戻ってきた後方散乱光は散乱回数が多くなり、光強度角度分布が等方的となり光伝播経路や光路長の情報は失われている。これに対して、ＴＯＦ情報には光伝播経路や光路長の情報が残っている。従って、深層領域の光学物性値分布の推定にはＴＯＦ情報を用い、浅層領域の光学物性値分布の推定には後方散乱光の角度分布情報を用いることが好ましい。なお、浅層領域及び深層領域からの後方散乱光の検出には、公知の様々な方法を用いてよい。
（実施例）
　図７は、実施例に係る光学物性値分布推定装置により計測を行う被計測物を表す図である。図７に示すように、実施例として、光学物性値分布推定装置１を用いて、第１層と第２層とで光学物性値が異なる被計測物２の光学物性値分布の計測を行う。第１層の厚さは２ｍｍ、被計測物２全体の厚さは７ｍｍ、幅は４０ｍｍである。
　図８は、図７に示す被計測物の光学物性値を表す図である。図８に示すように、第１層と第２層とにおいて、吸収係数μａ［１／ｍｍ］と散乱異方性ｇとは等しいが、散乱係数μｓ［１／ｍｍ］が異なる。なお、本実施例では、吸収がない（μａ＝０）散乱体を被計測物２として散乱係数分布を推定するが、吸収がある散乱体を被計測物２としても同様に散乱係数分布と吸収係数分布とを推定することが可能である。
　図９は、ＴＯＦ計測部による計測結果を表す図である。図９に示すように、照明部１１と受光部１２との間の距離ｌ１［ｍｍ］（図２参照）を大きくすると、照明部１１から照射する光と受光部１２が受光した光との強度変調の位相ずれφの量［ｐｓ］が大きくなる。ただし、位相ずれφは、その他の公知技術により計測してもよい。
　続いて、図９の計測値を用いた逆解析演算により、被計測物２の散乱係数分布を推定する。図１０は、図７に示す被計測物をメッシュ状に分割する様子を表す図である。図１０に示すように、被計測物２をメッシュ状に分割する。図１１は、図１０に示すメッシュの拡大図である。図１１に示すように、各メッシュは、さらに角度θ（２２．５°）で１６等分された領域に分割されている。
　図１２は、図１０に示す被計測物に照明部と受光部とを配置した様子を表す図である。図１２に示すように、被計測物２の中央に配置された照明部１１は、被計測物２の表面に直交するように照射角度９０°で光を照射する。そして、照明部１１から照射された光は、±１１．２５°で拡がるものとする。また、被計測物２の表面全体に配置された受光部１２は、被計測物２の表面に直交する受光角度９０°から±９０°の全方位からの光を受光する。
　この状態で、公知の光輸送方程式及びＰＣＯ（偏微分方程式を制約条件とした最適化）アルゴリズム（例えば、学術文献「Ｏｐｔｉｃａｌ　ｔｏｍｏｇｒａｐｈｙ　ａｓ　ａ　ＰＤＥ－ｃｏｎｓｔｒａｉｎｅｄ　ｏｐｔｉｍｉｚａｔｉｏｎ　ｐｒｏｂｌｅｍ」　Ｉｎｖｅｒｓｅ　Ｐｒｏｂｌｅｍｓ　２１　（２００５）　１５０７－１５３０を参照）により予測値を推定する。ただし、光学物性値分布の推定方法は特に限定されない。例えば、光輸送方程式に代えて拡散方程式を用い光学物性値分布を推定してもよい。また、例えば、光学物性値分布を仮定して光伝播モデルを用いて数値計算を事前に行い、光学物性値分布とこれに対応する後方散乱光の計算値とをルックアップテーブルとしておき、Ｓｉｍｕｌａｔｅｄ　Ａｎｎｅａｌｉｎｇ法等のアルゴリズムを用いて光学物性値分布を推定してもよい。
　図１３は、第１の推定工程による推定結果を表す図である。図１３に示すように、例えば、散乱係数の初期値を１に仮定して逆解析演算を行うことにより、散乱係数の推定結果が算出される。この推定結果は、深層（第２層）において、実際の値と略等しい。
　図１４は、第２の推定工程による初期値の設定方法を説明するための図である。図１４に示すように、照明部１１の直下の柱状の領域Ａの推定結果により領域Ａ以外の領域の光学物性値を置換して第２の工程の初期値とする。換言すると、第２の光学物性値分布は、第１の光学物性値分布のうち、被計測物２に照射した光に沿って被計測物２の表面に垂直な方向の垂直領域Ａの光学物性値により、被計測物２の全領域の光学物性値を置換した値を初期値として推定される。この領域Ａは、照明部１１の直下であり光学物性値が最も高精度に推定されている領域だからである。このように、第２の推定工程において、第１の光学物性値分布の一部を初期値とした逆解析演算を行ってもよい。また、この場合、第１の推定工程において、第２の推定工程に用いる光学物性値を優先的に算出し、第１の推定工程におけるその他の光学物性値の算出と、第２の推定工程とを並行して行ってもよい。なお、被計測物２の内部の構造によっては、照明部１１からの光を被計測物２の表面に対して斜めに照射した方がよい場合もある。この場合には、照明部１１から照射する光に沿った柱状の領域の推定結果によりこの領域以外の領域の光学物性値を置換して第２の工程の初期値とすればよい。また、照明部１１と受光部１２とが１個だけでなく、複数個の光源及び光検出器を有する場合もある。その場合には、複数個の光源から照射する光に沿った複数の柱状の領域の各推定結果に基づいて、これらの領域以外の領域の光学物性値を置換して第２の推定工程の初期値とすればよい。なお、第２の推定工程の初期値として、複数の柱状の領域の各推定結果のいずれか１つを用いてもよいし、各柱状の領域の推定結果の平均や重み付けをして算出した値等、公知の様々な方法により算出された値を用いてもよい。
　図１５～図１８は、角度分布計測部による計測結果を表す図である。図１５～図１８は、照明部２１と受光部２２との距離ｌ２（図２参照）を０ｍｍ～１．０ｍｍまで変化させ、各距離ｌ２において、受光部２２の受光角度を０°～６０°まで変化させた場合の計測値の一部を表すものである。このように、照明部２１と受光部２２との間の距離ｌ２と受光角度とを変化させると、光強度の角度分布が計測できる。なお、受光角度は、受光部２２が受光する光の中心角であり、受光部２２は、受光角度の±１１．２５°の光を受光するものとした。
　そして、図１５～図１８の計測値を用いて、図１２に示す計算モデルを用いることにより被計測物２内を光が伝播する様子を算出し、光学物性値を算出する。なお、ここでは、被計測物２の表面全体に配置された受光部２２は、計測値と対応するように受光角度を０°～６０°まで変化させて計算を行う。
　図１９は、第２の推定工程による推定結果を表す図である。図１９に示すように、第１の推定工程の値を初期値として光学物性値を推定することにより、実際の値と略等しい散乱係数の推定結果を得ることができた。
　なお、光強度角度分布を計測する際に、被計測物２の第１層の概略の等価散乱係数μｓ’ｒｏｕｇｈが文献等により既知である場合には、第１層に十分光が到達するように照射角度及び受光角度を設定すればよい。図２０は、第１層の厚さを推定する方法の一例を表す図である。図２０に示すように、被計測物２の第１層の概略の等価散乱係数μｓ’ｒｏｕｇｈが既知である場合には、第１層の厚さは例えば１／μｓ’ｒｏｕｇｈにより算出することができる。このように、被計測物２の表層領域における光学物性値が小さい領域の深さを予測し、その深さに応じて被計測物２に光を照射する角度、及び被計測物２からの後方散乱光を受光する角度を定め、定められた角度で第２の計測値を計測することによって、より短時間で高精度に光学物性値を推定することができる。すなわち、第２の計測値は、既知である散乱の強さを表す光学物性値（散乱係数等）に応じて定められた計測方法により計測されることが好ましい。
　また、上述した実施の形態において、ＴＯＦ計測部３２が、ＴＯＦ情報として、照明部１１から照射する光と受光部１２が受光した光との強度変調の位相ずれを算出する構成を説明したがこれに限らない。例えば、パルスレーザーである照明部１１からパルス光を照射し、ストリークカメラである受光部１２が被計測物２における後方散乱光の時間分解波形を検出する場合、ＴＯＦ計測部３２は、ＴＯＦ情報として、照明部から照射したパルス光に対する受光部１２が受光した光の時間遅れを算出してもよい。
　さらなる効果や変形例は、当業者によって容易に導き出すことができる。よって、本発明のより広範な態様は、以上のように表し、かつ記述した特定の詳細及び代表的な実施の形態に限定されるものではない。従って、添付のクレーム及びその均等物によって定義される総括的な発明の概念の精神又は範囲から逸脱することなく、様々な変更が可能である。
　１　光学物性値分布推定装置
　２　被計測物
　１０　ＴＯＦセンサユニット
　１１、２１　照明部
　１２、１２ａ、１２ｂ、２２、２２ａ、２２ｂ　受光部
　２０　角度分布計測装置
　２０Ａ　角度分布計測照明ユニット
　２０Ｂ　角度分布計測受光ユニット
　３０　処理装置
　３１　制御部
　３２　ＴＯＦ計測部
　３３　角度分布計測部
　３４　逆解析演算部
　３５　画像生成部
　３６　記憶部
　４０　入力装置
　５０　表示装置

Claims (14)

  1. 　被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置が行う光学物性値分布の推定方法であって、
    　前記被計測物に照射した光の等方的な後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光学物性値分布である第１の光学物性値分布を推定する第１の推定工程と、
    　前記被計測物に照射した光の前記第１の計測値より非等方的な後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の光学物性値分布である第２の光学物性値分布を推定する第２の推定工程と、
    　を含むことを特徴とする光学物性値分布の推定方法。
  2. 　第２の光学物性値分布は、前記被計測物内部の前記第１の光学物性値より浅層の光学物性値分布であることを特徴とする請求項１に記載の光学物性値分布の推定方法。
  3. 　前記第２の推定工程は、前記第１の光学物性値分布のうち、前記被計測物に照射した光に沿って前記被計測物の表面に垂直な方向の垂直領域の光学物性値により、前記被計測物の全領域の光学物性値を置換した値を初期値として前記第２の光学物性値分布を推定することを特徴とする請求項１又は２に記載の光学物性値分布の推定方法。
  4. 　前記第１の計測値は、前記被計測物にパルス光又は周期的に強度を変調させた光を照射して計測したＴＯＦ情報を含む値であり、
    　前記第１の推定工程は、前記ＴＯＦ情報に基づいて前記第１の光学物性値分布を推定することを特徴とする請求項１～３のいずれか１つに記載の光学物性値分布の推定方法。
  5. 　前記第２の計測値は、前記被計測物に光を照射する角度、又は前記被計測物からの後方散乱光を受光する角度の少なくとも一方を、前記被計測物の表面に対して鋭角に角度をつけて計測した値であり、
    　前記第２の推定工程は、前記第２の計測値の光強度角度分布に基づいて前記第２の光学物性値分布を推定することを特徴とする請求項１～４のいずれか１つに記載の光学物性値分布の推定方法。
  6. 　前記第１の計測値は、光音響現象、直交ニコル法の少なくともいずれか１つを用いて計測した値であることを特徴とする請求項１～５のいずれか１つに記載の光学物性値分布の推定方法。
  7. 　前記第２の計測値は、ＯＣＴ法、平行ニコル法の少なくともいずれか１つを用いて計測した値であることを特徴とする請求項１～６のいずれか１つに記載の光学物性値分布の推定方法。
  8. 　前記第２の計測値は、既知である散乱の強さを表す光学物性値に応じて定められた計測方法により計測した値であることを特徴とする請求項１～７のいずれか１つに記載の光学物性値分布の推定方法。
  9. 　前記計測方法は、前記被計測物の表層領域における前記光学物性値が小さい領域の深さを予測し、該深さに応じて前記被計測物に光を照射する角度、及び前記被計測物からの後方散乱光を受光する角度を定め、該定められた角度で前記第２の計測値を計測することを特徴とする請求項８に記載の光学物性値分布の推定方法。
  10. 　前記深さは、前記被計測物の概略の等価散乱係数μｓ’ｒｏｕｇｈを用いて、１／μｓ’ｒｏｕｇｈにより算出されることを特徴とする請求項９に記載の光学物性値分布の推定方法。
  11. 　被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置が行う光学物性値分布の推定方法であって、
    　前記被計測物に照射した光の後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光エネルギーの損失が大きいと推測される領域を含む第１の領域の光学物性値分布である第１の光学物性値分布を推定する第１の推定工程と、
    　前記被計測物に照射した光の後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の前記第１の領域より光エネルギーの損失が小さいと推測される第２の領域の光学物性値分布である第２の光学物性値分布を推定する第２の推定工程と、
    　を含むことを特徴とする光学物性値分布の推定方法。
  12. 　被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置が行う光学物性値分布の推定方法であって、
    　前記被計測物の第１の領域に照射した光の後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光学物性値分布である第１の光学物性値分布を推定する第１の推定工程と、
    　前記第１の領域に包含される第２の領域に照射した光の後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の前記第１の光学物性値分布より浅層の光学物性値分布である第２の光学物性値分布を推定する第２の推定工程と、
    　を含むことを特徴とする光学物性値分布の推定方法。
  13. 　被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置が行う光学物性値分布の推定プログラムであって、
    　前記被計測物に照射した光の等方的な後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光学物性値分布である第１の光学物性値分布を推定する第１の推定ステップと、
    　前記被計測物に照射した光の前記第１の計測値より非等方的な後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の前記第１の光学物性値分布より浅層の光学物性値分布である第２の光学物性値分布を推定する第２の推定ステップと、
    　をコンピュータに実行させることを特徴とする光学物性値分布の推定プログラム。
  14. 　被計測物の内部の光学物性値分布を推定する光学物性値分布推定装置であって、
    　前記被計測物に照射した光の等方的な後方散乱光を計測した第１の計測値を記憶部から読み出し、逆解析演算により前記被計測物内部の光学物性値分布である第１の光学物性値分布を推定し、
    　前記被計測物に照射した光の前記第１の計測値より非等方的な後方散乱光を計測した第２の計測値を記憶部から読み出し、前記第１の光学物性値分布の少なくとも一部を初期値とした逆解析演算により、前記被計測物内部の前記第１の光学物性値分布より浅層の光学物性値分布である第２の光学物性値分布を推定する逆解析演算部を備えることを特徴とする光学物性値分布推定装置。
PCT/JP2017/014869 2017-04-11 2017-04-11 光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置 WO2018189815A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/JP2017/014869 WO2018189815A1 (ja) 2017-04-11 2017-04-11 光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置
US16/594,166 US10955337B2 (en) 2017-04-11 2019-10-07 Method of estimating optical physical property value distribution, computer-readable recording medium and optical physical property value distribution estimating apparatus

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2017/014869 WO2018189815A1 (ja) 2017-04-11 2017-04-11 光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US16/594,166 Continuation US10955337B2 (en) 2017-04-11 2019-10-07 Method of estimating optical physical property value distribution, computer-readable recording medium and optical physical property value distribution estimating apparatus

Publications (1)

Publication Number Publication Date
WO2018189815A1 true WO2018189815A1 (ja) 2018-10-18

Family

ID=63792729

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/014869 WO2018189815A1 (ja) 2017-04-11 2017-04-11 光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置

Country Status (2)

Country Link
US (1) US10955337B2 (ja)
WO (1) WO2018189815A1 (ja)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2013019696A (ja) * 2011-07-07 2013-01-31 Univ Of Tsukuba 生体計測装置および画像作成方法
JP2013101060A (ja) * 2011-11-09 2013-05-23 Canon Inc 被検体光計測装置及び被検体光計測方法
JP2014027960A (ja) * 2012-06-29 2014-02-13 Univ Of Tsukuba 生体計測装置の計測データ選択方法、生体計測装置の光出射位置決定方法、および生体計測装置

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6662030B2 (en) * 1998-05-18 2003-12-09 Abbott Laboratories Non-invasive sensor having controllable temperature feature
NO325061B1 (no) 2001-03-06 2008-01-28 Photosense As Fremgangsmate og arrangement for bestemmelse av den optiske egenskap av et multisjiktvev
DE10163972B4 (de) 2001-12-22 2005-10-27 Roche Diagnostics Gmbh Verfahren und Vorrichtung zur Bestimmung eines Lichttransportparameters und eines Analyten in einer biologischen Matrix
JP4327738B2 (ja) 2005-01-18 2009-09-09 株式会社東芝 生体光計測装置及び生体光計測方法
JP4664260B2 (ja) * 2005-09-21 2011-04-06 シャープ株式会社 表示装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2013019696A (ja) * 2011-07-07 2013-01-31 Univ Of Tsukuba 生体計測装置および画像作成方法
JP2013101060A (ja) * 2011-11-09 2013-05-23 Canon Inc 被検体光計測装置及び被検体光計測方法
JP2014027960A (ja) * 2012-06-29 2014-02-13 Univ Of Tsukuba 生体計測装置の計測データ選択方法、生体計測装置の光出射位置決定方法、および生体計測装置

Also Published As

Publication number Publication date
US20200033262A1 (en) 2020-01-30
US10955337B2 (en) 2021-03-23

Similar Documents

Publication Publication Date Title
JP5463545B2 (ja) 濃度定量装置、濃度定量方法及びプログラム
JP5782314B2 (ja) 生体計測装置および画像作成方法
US20170071476A1 (en) Image generating apparatus, image generating method, and program
US9538970B2 (en) Generating attenuation image data and phase image data in an X-ray system
US9833179B2 (en) Blood component analyzing method and blood component analyzing apparatus
EP2302362B1 (fr) Dispositif et procédé de reconstruction spatiale d'une cartographie de fluorescence
CN103300880A (zh) 被检体信息获得装置和被检体信息获得方法
JP2018532121A (ja) 波長可変レーザを用いた大面積ｏｃｔシステムと３次元イメージ補正方法
JP6362420B2 (ja) 被検体情報取得装置、被検体情報取得方法、および、プログラム
RU2602718C2 (ru) Устройство для получения информации о субъекте и способ получения информации в отношении субъекта
US10034611B2 (en) Subject information obtaining apparatus and subject information obtaining method
JP2017026611A (ja) 光検出装置および光検出方法
JP5183406B2 (ja) 生体情報処理装置及び生体情報処理方法
JP6358573B2 (ja) 乳房計測装置の作動方法及び乳房計測装置
JP6894088B2 (ja) 散乱体濃度計測装置及びその方法
TW201935035A (zh) 距離飛行時間模組
JP2013070822A (ja) 濃度定量装置、光吸収係数算出方法、等価散乱係数算出方法、濃度定量方法、光吸収係数の算出を行うプログラム及び濃度の算出を行うプログラム
JP2019201760A (ja) 血管検知装置及びその方法
WO2018189815A1 (ja) 光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置
JP2008155011A (ja) 密度計測装置およびその方法
JP2008527321A (ja) 混濁媒質の光学特性を決定するための方法
JP2014073300A (ja) 被検体情報取得装置および被検体情報取得装置の制御方法
US10921248B2 (en) Measurement apparatus and measurement method
US20220361806A1 (en) Method For Detecting Demineralization Of Tooth Substance
WO2020246455A1 (ja) 無侵襲計測装置、方法およびプログラム

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

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 17905438

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: JP