WO2013005635A1 - 生体計測装置および画像作成方法 - Google Patents

生体計測装置および画像作成方法 Download PDF

Info

Publication number
WO2013005635A1
WO2013005635A1 PCT/JP2012/066552 JP2012066552W WO2013005635A1 WO 2013005635 A1 WO2013005635 A1 WO 2013005635A1 JP 2012066552 W JP2012066552 W JP 2012066552W WO 2013005635 A1 WO2013005635 A1 WO 2013005635A1
Authority
WO
WIPO (PCT)
Prior art keywords
pixel
reconstructed image
image
measurement site
light
Prior art date
Application number
PCT/JP2012/066552
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 EP12807701.3A priority Critical patent/EP2730231B1/en
Priority to US14/130,690 priority patent/US9282932B2/en
Priority to CN201280033086.7A priority patent/CN103635148B/zh
Publication of WO2013005635A1 publication Critical patent/WO2013005635A1/ja

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • 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
    • A61B10/0041—Detection of breast cancer
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0073—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence by tomography, i.e. reconstruction of 3D images from 2D projections
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0082—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence adapted for particular medical purposes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • 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/5229—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data combining image data of a patient, e.g. combining a functional image with an anatomical image
    • A61B6/5247—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data combining image data of a patient, e.g. combining a functional image with an anatomical image combining images from an ionising-radiation diagnostic technique and a non-ionising radiation diagnostic technique, e.g. X-ray and ultrasound
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/13—Tomography
    • A61B8/15—Transmission-tomography
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
    • 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

Definitions

  • the present invention relates to a biological measurement apparatus and an image creation method.
  • a so-called diffuse optical tomography is a device that noninvasively measures internal information of a living body such as the head and breast, and obtains internal information using the light absorption characteristics of the living body.
  • Tomography has been proposed.
  • light is irradiated from a predetermined irradiation position to a part of a living body to be measured, and light propagated while being scattered inside the part is detected at a predetermined detection position. From the measurement results such as intensity and time waveform, it is possible to obtain internal information of the part, that is, information on a light absorber such as a tumor existing inside the part.
  • Patent Document 1 describes a living body measurement method using diffuse optical tomography.
  • Non-Patent Documents 1 and 2 describe a successive approximation image reconstruction method in diffuse optical tomography.
  • FIG. 17 is a schematic diagram for explaining such a phenomenon, and shows a measurement site 100 and a light irradiation unit 101 and a light detection unit 102 provided on the surface of the measurement site 100.
  • the shorter the flight time of photons emitted from the light irradiation position 101 and reaching the light detection position 102 the shorter the flight distance, and the flight path is limited.
  • the longer the photon flight time the longer the flight distance and the flight path is not limited.
  • the data having a short flight time of the photon includes many flight paths R1 passing through the area A1 shown in FIG. 17, that is, the area close to the surface of the measurement target part. Further, the data having a long flight time of the photon includes many flight paths R2 that pass through the region A2 shown in FIG. 17, that is, the region far from the surface of the measurement target part. Therefore, the information amount in the region far from the surface of the measurement site is smaller than the information amount in the region near the surface of the measurement site. As a result, the spatial resolution and noise in the region far from the surface of the measurement site are larger than the spatial resolution and noise in the region near the surface of the measurement site.
  • the present invention has been made in view of such problems, and a living body measurement apparatus capable of suppressing the difference in spatial resolution and noise characteristics depending on the position inside the measurement target part and creating a more uniform image. It is another object of the present invention to provide an image creation method.
  • a first biological measurement apparatus includes a light irradiation unit that irradiates light to a measurement site of a subject, and light detection that detects diffused light from the measurement site. And a calculation unit that calculates a light absorption coefficient distribution in the measurement site based on an output signal from the light detection unit and creates a reconstructed image related to the inside of the measurement site.
  • a second living body measurement apparatus includes an irradiation unit that irradiates a measurement site of a subject with radiation or sound waves, a detection unit that detects diffused radiation or sound waves from the measurement site, and a detection A calculation unit that calculates a radiation or sound wave absorption coefficient distribution in the measurement site based on an output signal from the unit and creates a reconstructed image related to the inside of the measurement site.
  • J coefficients w j greater than 0 and less than or equal to 1 where subscript j is an integer from 1 to J, and J is the number of pixels in the reconstructed image).
  • a reconstructed image is created by performing a successive approximation operation using the above iterative formula (1).
  • radiation refers to an electromagnetic wave having a short wavelength such as X-ray, ⁇ -ray, or microwave
  • sound wave refers to a wave such as an ultrasonic wave.
  • successive approximation calculation for image reconstruction is performed using J coefficients w 1 to w J set for each pixel of the reconstructed image. For example, by setting these coefficients w 1 to w J so that the convergence speed of the pixel coincides with a region where the convergence speed is the slowest in N iterations, the convergence speed is made uniform, and the measurement target part It is possible to create a more uniform image while suppressing differences in spatial resolution and noise characteristics depending on the internal position.
  • the calculation unit has a convergence speed (hereinafter referred to as a convergence speed) of a partial area having the slowest convergence speed in N iterations among a plurality of partial areas each including a plurality of pixels formed by dividing a reconstructed image.
  • a convergence speed a convergence speed of a partial area having the slowest convergence speed in N iterations among a plurality of partial areas each including a plurality of pixels formed by dividing a reconstructed image.
  • C N is calculated, and M numerical values v m satisfying 0 ⁇ v m ⁇ 1 are prepared, where m is an integer from 1 to M, and for each m from 1 to M, Next iteration Calculating a pixel value x 1 of the J each pixel of the (N) ⁇ x J (N ) by performing N iterations operation using, from the pixel value x 1 (N) ⁇ x J (N )
  • a numerical value v m when the convergence speed of each partial area obtained substantially coincides with the minimum convergence speed C N may be set as the coefficient w j of a plurality of pixels included in the partial area. Thereby, the effect mentioned above can be acquired more suitably.
  • the biological measurement apparatus may be configured such that the calculation unit sets the coefficient w j of the partial region with the slowest convergence speed to 1 in N iterations. Thereby, the effect mentioned above can be acquired more suitably.
  • a first image creation method irradiates a measurement site of a subject with light, detects diffused light from the measurement site, and based on the detection signal, a light absorption coefficient in the measurement site
  • Subscript j is an integer from 1 to J, where J is the number of pixels of the reconstructed image) (Where k is an integer from 1 to N, where N is the number of iterations.
  • X j (k) is the pixel value at the k-th iteration of the j-th pixel
  • d j ( k) is an update amount at the time of the k-th iterative calculation of the j-th pixel.
  • the second image creation method irradiates the measurement site of the subject with radiation or sound waves, detects diffused radiation or sound waves from the measurement site, and detects the subject based on the detection signal.
  • a reconstructed image is created by performing an operation.
  • radiation refers to an electromagnetic wave having a short wavelength such as X-ray, ⁇ -ray, or microwave
  • sound wave refers to a wave such as an ultrasonic wave.
  • successive approximation calculation for image reconstruction is performed using J coefficients w 1 to w J set for each pixel of the reconstructed image. For example, by setting these coefficients w 1 to w J so that the convergence speed of the pixel coincides with a region where the convergence speed is the slowest in N iterations, the convergence speed is made uniform, and the measurement target part It is possible to create a more uniform image while suppressing differences in spatial resolution and noise characteristics depending on the internal position.
  • the image creation method includes a convergence speed (hereinafter referred to as a minimum convergence speed) of a partial area having the slowest convergence speed in N iterations among a plurality of partial areas each formed by dividing a reconstructed image and including a plurality of pixels.
  • a convergence speed hereinafter referred to as a minimum convergence speed
  • C N is determined, and M numerical values v m satisfying 0 ⁇ v m ⁇ 1 are prepared (where m is an integer from 1 to M), and for each m from 1 to M, the following iterative formula Calculating a pixel value x 1 of the J each pixel of the (N) ⁇ x J (N ) by performing N iterations operation using, from the pixel value x 1 (N) ⁇ x J (N )
  • a numerical value v m when the convergence speed of each partial area obtained substantially coincides with the minimum convergence speed C N may be set as the coefficient w j of a plurality of pixels included in the partial area. Thereby, the effect mentioned above can be acquired more suitably.
  • the image creation method may be configured such that the coefficient w j of the partial region having the slowest convergence speed in N iterations is set to 1. Thereby, the effect mentioned above can be acquired more suitably.
  • the living body measurement apparatus and the image creation method of the present invention it is possible to create a more uniform image by suppressing the difference in spatial resolution and noise characteristics depending on the position inside the measurement site.
  • FIG. 1 is a diagram illustrating a configuration of a biological measurement apparatus according to an embodiment.
  • FIG. 2 is a flowchart showing a specific method for determining a coefficient (step size).
  • FIG. 3 is a diagram illustrating two types of measurement sites to be reconfigured in the simulation.
  • FIG. 4 is a chart showing conditions in the forward problem analysis for simulating the flight of photons inside the measurement site.
  • FIG. 5 is a chart showing conditions in inverse problem analysis for reconstructing an image from detected photon histograms.
  • FIG. 6 shows four images used to determine the step size in the simulation.
  • FIG. 7 shows (a) the image after reconstruction by simulation corresponding to the image shown in FIG. 3 (a), and (b) the step size corresponding to the image shown in FIG. 3 (a).
  • FIG. 8 is a graph showing changes in pixel values on the three lines shown in FIG.
  • FIG. 9 is a diagram illustrating three lines assumed on the image.
  • FIG. 10 shows step sizes corresponding to (a) the image after reconstruction by simulation corresponding to the image shown in FIG. 3 (b), and (b) the image shown in FIG. 3 (b). It is the image after the reconstruction by the conventional method which is not used.
  • FIG. 11 is a graph showing changes in pixel values on the two lines shown in FIG.
  • FIG. 12 is a diagram showing two lines assumed on the image.
  • FIG. 13 shows (a) an image after reconstruction according to this embodiment when statistical noise is added, corresponding to the image shown in FIG. 3 (a), and (b) shown in FIG.
  • FIG. 3 (a) It is the image after the reconstruction by the conventional method at the time of adding statistical noise corresponding to the obtained image.
  • FIG. 14 shows (a) an image after reconstruction according to this embodiment when statistical noise is added, corresponding to the image shown in FIG. 3 (b), and (b) shown in FIG. 3 (b). It is the image after the reconstruction by the conventional method at the time of adding statistical noise corresponding to the obtained image.
  • FIG. 15 is a graph showing the NSD for each partial region.
  • FIG. 16 is a chart showing calculation results of standard deviation.
  • FIG. 17 is a schematic diagram for explaining a phenomenon in which a nonuniform image is generated in diffuse optical tomography.
  • FIG. 1 is a diagram showing a configuration of a biological measurement apparatus 10 according to an embodiment of the present invention.
  • the living body measurement apparatus 10 of the present embodiment irradiates light to a measurement site B of a subject who is a measurement target, detects diffused light (return light), and detects the detected position and light amount data (for example, time)
  • It is a device that estimates the average flight path and average optical path length of photons based on the decomposition photon histogram) and images in-vivo information as an image reconstruction problem.
  • the image obtained by this apparatus is a functional image of a body tissue, for example, which visualizes the position of a tumor and the distribution of oxygenated hemoglobin and deoxygenated hemoglobin.
  • a to-be-measured part B a head, a female breast, etc. are assumed, for example.
  • the biological measurement apparatus 10 includes a light irradiating unit that irradiates measurement light into the measurement site B, a light detection unit that detects diffused light generated from the measurement site B due to light irradiation from the light irradiation unit, and a light A calculation unit that calculates a spatial distribution of the absorption coefficient of the measurement site B based on an output signal from the detection unit and creates a reconstructed image of the measurement site B;
  • the light irradiating unit of the present embodiment includes a light emitting end, a light source 22, and an optical switch 24 that each of the n light emitting / detecting ends 16 attached to the measurement site B have.
  • the light source 22 for example, a laser diode can be used.
  • the wavelength of the measurement light is preferably a wavelength in the near infrared region of about 700 nm to 900 nm from the relationship between the transmittance of the living body and the absorption coefficient of the absorber to be quantified.
  • Measured light is emitted from the light source 22 as continuous light, for example.
  • the measurement light emitted from the light source 22 is irradiated to the measurement site B from the light emission / detection end 16.
  • the optical switch 24 is a 1-input n-output optical switch, which inputs light from the light source 22 through the light source optical fiber 26 and supplies the light to each of the n light emitting / detecting ends 16 in order. To do. That is, the optical switch 24 sequentially selects the n output optical fibers 28 connected to each light output / detection end 16 one by one, and optically connects the output optical fibers 28 and the light source 22. .
  • the light detection unit of the present embodiment includes a light detection end included in each of the n light emission / detection ends 16 described above, n light detectors 30 corresponding to each of the n light emission / detection ends 16, and The n shutters 32 are arranged in front of the input unit of each photodetector. Diffused light from the measurement site B incident on the light detection end of each light emission / detection end 16 is input to each of the n light detectors 30 via the detection optical fiber 34. The photodetector 30 generates an analog signal according to the light intensity of the diffused light that has reached the corresponding light emitting / detecting end 16.
  • the photodetector 30 various things such as a photomultiplier tube (PMT: Photomultiplier Tube), a photodiode, an avalanche photodiode, and a PIN photodiode can be used.
  • PMT Photomultiplier Tube
  • a photodiode an avalanche photodiode
  • a PIN photodiode a photodiode
  • a signal processing circuit 36 is connected to the signal output terminal of the photodetector 30, and the signal processing circuit 36 A / D converts the analog signal output from the photodetector 30 in accordance with the light intensity of the diffused light.
  • the digital signal is generated and provided to the arithmetic unit 14.
  • the calculation unit 14 calculates the light absorption coefficient distribution in the measurement site B and creates a reconstructed image relating to the inside of the measurement site B.
  • the calculation unit 14 is, for example, a CPU (Central Processing). Unit) and a computer having storage means such as a memory.
  • the calculation unit 14 may further have a function of controlling light emission of the light source 22, operation of the optical switch 24, and opening / closing of the shutter 32.
  • a recording / display unit 38 is connected to the calculation unit 14, and a calculation result in the calculation unit 14, that is, a reconstructed image of the measurement site B can be visualized.
  • Calculation of internal information of the measurement site B is performed as follows, for example.
  • the measurement light is sequentially irradiated from each of the n light emission / detection ends 16 to the inside of the measurement site B, and the light that has passed through the measurement site B and diffused is passed through the n light emission / detection ends 16. And detected by n photodetectors 30. Based on the detection result, the spatial distribution of the absorption coefficient inside the measurement site B is calculated, and a reconstructed image including information (internal information) on the position and shape of the absorber such as a tumor is created.
  • image reconstruction is performed by solving the maximization problem of the above equation (6) by the gradient method using the OS-Convex method.
  • the OS-Convex method measurement data is divided into partial data sets called subsets, and the solution is updated for each subset so that the evaluation function corresponding to each subset increases.
  • the solution can be updated by the number of subsets in one iterative operation, so that the convergence speed is improved.
  • k is an integer from 1 to N (N is the number of iterations), the number of subsets is Q, q is an integer from 1 to Q, and the data set of the q-th subset is defined as Sq.
  • a specific iteration formula of the OS-Convex method is expressed by the following equations (8) and (9).
  • an iterative formula of a general successive approximation method is represented by the following formula (10).
  • the (k + 1) -th value x j (k + 1) is represented by the sum of the k-th value x j (k) and the update amount d j (k) .
  • the update amount d j (k) is expressed by the following equation (11).
  • the spatial resolution (image quality) of the reconstructed image after N iterations can be made uniform by bringing the convergence speeds of the values x 1 (N) to x J (N) in J pixels close to each other. it can.
  • the update amount d j (k) is the maximum update amount that increases the evaluation function in the j-th pixel, it is larger than 1 with respect to the update amount d j (k) of the pixel having a slow convergence speed.
  • Multiplying by a coefficient is not preferable because the evaluation function does not increase but decreases. Therefore, by multiplying the update amount d j (k) of a pixel having a high convergence speed by a coefficient smaller than 1, it is possible to increase the evaluation function and to adjust the convergence speed to a pixel having a low convergence speed. It becomes.
  • a coefficient w j (hereinafter referred to as a step size) satisfying the following equation (12) is set for each pixel, and the step size w j is multiplied by the update amount d j (k) , whereby each pixel value x It becomes possible to arbitrarily control the convergence speed of 1 (N) to x J (N) .
  • the step size w j it is necessary to determine the pixel having the slowest convergence speed after obtaining the convergence speed for each pixel based on the measurement data.
  • the convergence speed for each pixel is the contrast recovery rate (CRC: Contrast Recovery). Coefficient).
  • an image is divided into a plurality of partial areas each including a plurality of pixels, a spot area is arranged in each partial area, and the CRC at the spot is set as a value corresponding to the convergence speed of the area.
  • the CRC is defined by the following equation (15). However, SP represents a pixel value in the spot area, and BG represents a pixel value outside the spot area (background area). Further, the subscript m indicates an average value in the region, the subscript R indicates a pixel value after image reconstruction, and the subscript Tr indicates a pixel value based on measurement data.
  • FIG. 2 is a flowchart showing a specific method for determining the step size w j based on the convergence speed evaluation method described above.
  • this method first, an image is divided into E partial areas (step S1).
  • a set of pixel values in a certain arbitrary partial area e is denoted as R e .
  • Step S2 an arbitrary number of iterations N is set, a partial region having the smallest CRC (that is, the slowest convergence speed) is determined in N iterations, and the CRC is defined as the lowest convergence speed C N (Ste S2).
  • M arbitrary numerical values v m satisfying 0 ⁇ v m ⁇ 1 are prepared (where m is an integer from 1 to M, and v 1 to v M are different from each other) (step S3). ). Note that step S3 may be performed before step S2 or step S1.
  • step S4 the pixel values x 1 (N) to x J (N ) of the J pixels are obtained by performing N iterations for each m from 1 to M using the following iterative formula (16). ) And a reconstructed image is created (step S4).
  • the convergence speed of a certain partial area obtained from the pixel values x 1 (N) to x J (N) substantially coincides with the minimum convergence speed C N the numerical value v m at that time is stored in the partial area.
  • a step size w j of a plurality of pixels included is set (step S5). Thereafter, steps S4 and S5 described above are repeated M times, or steps S4 and S5 are repeated until the step size w j is determined for all partial regions except the partial region having the slowest convergence speed.
  • step S6 the step size w j of the partial region having the slowest convergence speed in N iterations is set to 1 (step S6).
  • the step size w j is determined for all the pixels in all the partial areas.
  • Diffuse light tomography is an image reconstruction method using a diffusive flight path of near-infrared light in a living body.
  • Diffuse optical tomography image reconstruction uses the successive approximation image reconstruction method rather than the analytical method, but the normal successive approximation image reconstruction method such as the conjugate gradient method uses spatial resolution and noise in the reconstruction image. There is a problem that non-uniform images with large variations in characteristics are generated.
  • the cause of the non-uniformity in image quality is that in diffuse optical tomography, the amount of information contained in measurement data varies greatly depending on the position inside the measurement site. That is, in diffuse optical tomography using the time-resolved measurement method, photons incident from the light incident end on the surface of the measurement site fly while repeatedly scattering within the measurement site and reach the light detection end of the surface of the measurement site. Is detected. The shorter the flight time of the photon from this incidence to detection, the shorter the flight distance and the more limited the flight path. Conversely, the longer the flight time, the longer the flight distance and the flight path is not limited. For this reason, the amount of information included in the measurement data changes according to the flight time of photons.
  • the ratio of the flight path close to the surface of the measurement site increases, and when the flight time is long, it is far from the surface of the measurement site (of the measurement site). Increases the percentage of flight paths (near the center). For this reason, the amount of information near the center of the measurement site is smaller than that near the surface. This causes variations in spatial resolution and noise in the reconstructed image such that the spatial resolution and noise in the peripheral portion of the reconstructed image are larger than the spatial resolution and noise in the vicinity of the center of the reconstructed image.
  • successive approximation calculation for image reconstruction is performed using J coefficients w 1 to w J set for each pixel of the reconstructed image. Then, as in this embodiment, these coefficients w 1 to w J are set so that the convergence speed of the pixel matches the partial area where the convergence speed is the slowest in N iterations. Can be made uniform, and it becomes possible to create a more uniform image by suppressing the difference in spatial resolution and noise characteristics depending on the position inside the measurement site.
  • FIG. 3A and FIG. 3B are diagrams showing two types of measurement sites to be reconstructed in this simulation, and each figure shows a plurality of light absorbers D such as tumors, for example.
  • the center coordinates (x, y) of the three light absorbers D shown in FIG. 3A are (21, 20), (64, 106), and (75, 51), respectively.
  • the center coordinates (x, y) of the four light absorbers D shown in FIG. 3B are (56, 79), (66, 24), (109, 49), and (109, 110), respectively. ).
  • the conditions common to the above steps (1) to (3) are as follows. First, the conditions in two types of analysis, forward problem analysis (Forward problem) that simulates the flight of photons inside the measurement site, and inverse problem analysis (Inverse problem) that reconstructs an image from detected photon histograms, These are shown in FIGS. 4 and 5, respectively. As shown in FIGS. 4 and 5, the grid size and the image size are different between the forward problem analysis and the inverse problem analysis because the continuous system is converted to the discrete system. Photons detected by the light detection unit are always acquired as discrete data. On the other hand, the value of the numerical simulation is calculated in a form close to a continuous function. Therefore, in this simulation, the continuous system is converted to the discrete system by downsampling.
  • Forward problem forward problem
  • Inverse problem inverse problem analysis
  • An image including pixels of 132 rows and 132 columns based on measurement data was divided into 12 rows and 12 columns of partial areas, and a hot spot was arranged at the center of each partial area.
  • the CRC of each partial region is obtained by reconstructing the four images shown in FIGS. 6A to 6D, and the step sizes w 1 to w J of all the pixels are obtained by the method described above. It was determined.
  • step size image having step sizes w 1 to w J as pixel values is an average value.
  • the filter was used to blur and obtain a step size image that changed smoothly.
  • an average value filter having a kernel size of 9 ⁇ 9 was used.
  • FIG. 7A is an image after reconstruction by this simulation corresponding to the image shown in FIG.
  • FIG. 7B is an image after reconstruction by a conventional method that does not use step sizes w 1 to w J
  • FIG. 8A to FIG. 8C are graphs showing changes in pixel values on the three lines L1 to L3 passing through the light absorber D as shown in FIG. A change in pixel value by simulation is shown, and a graph G22 shows a change in pixel value by a conventional method.
  • FIG. 10A is an image after reconstruction by this simulation corresponding to the image shown in FIG.
  • FIG. 10B is an image after reconstruction by a conventional method that does not use step sizes w 1 to w J corresponding to the image shown in FIG.
  • FIG. 10A is an image after reconstruction by this simulation corresponding to the image shown in FIG.
  • FIG. 10B is an image after reconstruction by a conventional method that does not use step sizes w 1 to w J corresponding to the image shown in FIG.
  • FIG. 11A and FIG. 11B are graphs showing changes in pixel values on two lines L4 and L5 passing through the light absorber D as shown in FIG. A change in pixel value by simulation is shown, and a graph G22 shows a change in pixel value by a conventional method.
  • the pixel value of the light absorber D in the peripheral portion of the image such as the upper left is relatively high, which is larger than the pixel value of the light absorber D near the center of the image. You can see that they are different.
  • the pixel values of the respective light absorbers D are the peripheral part and the center of the image. The value is almost uniform in the vicinity.
  • the change in pixel value in the peripheral portion of the image is large.
  • the change in the pixel value in the peripheral portion of the image is small.
  • the step size w j is increased near the center of the image and the update amount in the iterative calculation is increased, and the step size w j is decreased and the update amount is decreased in the peripheral portion of the image.
  • FIG. 13A is an image after reconstruction according to the present embodiment when statistical noise is added
  • corresponding to the image shown in FIG. 13B is an image after reconstruction by the conventional method when statistical noise is added
  • FIG. 14A shows an image after reconstruction according to the present embodiment when statistical noise is added
  • FIG. 14B is an image after reconstruction by the conventional method when statistical noise is added, corresponding to the image shown in FIG.
  • the conventional method that does not use the step sizes w 1 to w J increases the spatial resolution of the peripheral portion of the image and affects the noise.
  • the spatial resolution is not high and the influence of noise is small.
  • the update amount of the image peripheral portion in the iterative calculation is suppressed. Therefore, the influence of noise is remarkably suppressed, and the shape of the light absorber D existing at the upper left can be clearly recognized.
  • the influence of noise is suppressed to the same extent as in the conventional method.
  • the statistical noise near the center of the image is larger than the statistical noise at the periphery of the image because the resolution improvement speed varies depending on the presence or absence of the statistical noise, and it was determined on the assumption that there was no statistical noise. This is due to the use of size.
  • the biometric measurement method of the present embodiment has a change in pixel value as compared with the conventional method, but there is no significant change in the morphological characteristics of the statistical noise of the entire image. I understand that. That is, according to the living body measurement method of the present embodiment, it can be seen that the resolution improvement can be selectively delayed at the peripheral portion of the image.
  • FIG. 15 is a graph showing the NSD for each partial area, where the horizontal axis indicates the number of the partial area and the vertical axis indicates the value of NSD.
  • the subscript SD indicates a standard deviation.
  • a standard deviation using NSD for each partial area as a population is calculated, and this standard deviation is used as a value for evaluating noise non-uniformity between the partial areas.
  • FIG. 16 is a chart showing calculation results of standard deviation.
  • the living body measurement apparatus and the image creation method according to the present invention are not limited to the above-described embodiments, and various other modifications are possible.
  • the measurement site of the subject is irradiated with light, diffused light from the measurement site is detected, the light absorption coefficient distribution in the measurement site is calculated based on the detection signal, A reconstructed image relating to the inside of the measurement site is created.
  • the present invention is not limited to a living body measuring apparatus and an image creating method using light, but can also be applied to a living body measuring apparatus and an image creating method using radiation and sound waves.
  • the biological measurement apparatus includes an irradiation unit that irradiates a measurement site of a subject with radiation or sound waves, a detection unit that detects diffused radiation or sound waves from the measurement site, and a detection unit.
  • a calculation unit that calculates the absorption coefficient distribution of radiation or sound waves in the measurement site based on the output signal and creates a reconstructed image relating to the inside of the measurement site may be provided. In this case, the calculation unit may create a reconstructed image by the same method as in the above-described embodiment.
  • the image creating method irradiates a measurement site of a subject with radiation or sound waves, detects diffused radiation or sound waves from the measurement site, and based on the detection signal, A method may be used in which the absorption coefficient distribution of radiation or sound waves is calculated and a reconstructed image relating to the inside of the measurement site is created by the same method as in the above-described embodiment.
  • the radiation refers to an electromagnetic wave having a short wavelength such as X-ray, ⁇ -ray, or microwave
  • the sound wave refers to a wave such as an ultrasonic wave.
  • a light irradiation unit that irradiates light to a measurement site of a subject, a light detection unit that detects diffused light from the measurement site, and an output from the light detection unit
  • J coefficients w j greater than 0 and less than or equal to 1 where the subscript j is an integer from 1 to J, where J is the number of pixels in the reconstructed image
  • the reconstructed image is created by performing the successive approximation calculation using (1).
  • an irradiation unit that irradiates the measurement site of the subject with radiation or sound waves
  • a detection unit that detects diffused radiation or sound waves from the measurement site
  • detection A calculation unit that calculates a radiation or sound wave absorption coefficient distribution in the measurement site based on an output signal from the unit and creates a reconstructed image related to the inside of the measurement site.
  • J coefficients w j greater than 0 and less than or equal to 1 where subscript j is an integer from 1 to J, and J is the number of pixels in the reconstructed image).
  • a reconstructed image is created by performing successive approximation calculation using the above iterative formula (1).
  • radiation refers to an electromagnetic wave having a short wavelength such as X-rays, ⁇ -rays, or microwaves
  • sound waves refer to waves such as ultrasonic waves.
  • the calculation unit has a convergence speed (hereinafter referred to as a convergence speed) of a partial area having the slowest convergence speed in N iterations among a plurality of partial areas each formed by dividing a reconstructed image and including a plurality of pixels.
  • a convergence speed a convergence speed of a partial area having the slowest convergence speed in N iterations among a plurality of partial areas each formed by dividing a reconstructed image and including a plurality of pixels.
  • C N is calculated, and M numerical values v m satisfying 0 ⁇ v m ⁇ 1 are prepared, where m is an integer from 1 to M, and for each m from 1 to M,
  • the pixel values x 1 (N) to x J (N) of the J pixels are calculated by performing N iterations using the above iterative formula (2), and the pixel values x 1 (N) to A numerical value v m when the convergence speed of each partial area obtained from x J (N) substantially coincides with the minimum convergence speed C N may be set as a coefficient w j of a plurality of pixels included in the partial area. .
  • the biological measurement apparatus may be configured such that the calculation unit sets the coefficient w j of the partial region with the slowest convergence speed to 1 in N iterations. Thereby, the effect mentioned above can be acquired more suitably.
  • the measurement site of the subject is irradiated with light, diffused light from the measurement site is detected, and the light absorption coefficient in the measurement site is based on the detection signal.
  • the subscript j is an integer from 1 to J, where J is the number of pixels of the reconstructed image), and a reconstructed image is created by performing successive approximation using the above iterative formula (3) It is configured to do.
  • the measurement site of the subject is irradiated with radiation or sound waves, diffused radiation or sound waves from the measurement site are detected, and the target is detected based on the detection signal.
  • a reconstructed image is created by performing an operation.
  • radiation refers to an electromagnetic wave having a short wavelength such as X-rays, ⁇ -rays, or microwaves
  • sound waves refer to waves such as ultrasonic waves.
  • the image creation method includes a convergence speed (hereinafter referred to as a minimum convergence speed) of a partial area having the slowest convergence speed in N iterations among a plurality of partial areas each formed by dividing a reconstructed image and including a plurality of pixels.
  • a convergence speed hereinafter referred to as a minimum convergence speed
  • C N is calculated and M numerical values v m satisfying 0 ⁇ v m ⁇ 1 (where m is an integer from 1 to M) are prepared, and for each m from 1 to M, the above iterative formula ( 4), the pixel values x 1 (N) to x J (N) of the J pixels are calculated by performing N iterations, and the pixel values x 1 (N) to x J (N).
  • the numerical value v m when the convergence speed of each partial area obtained from (1 ) substantially coincides with the minimum convergence speed C N may be used as the coefficient w j of a plurality of pixels included in the partial area. Thereby, the effect mentioned above can be acquired more suitably.
  • the image creation method may be configured such that the coefficient w j of the partial region having the slowest convergence speed in N iterations is set to 1. Thereby, the effect mentioned above can be acquired more suitably.
  • the present invention can be used as a biological measurement apparatus and an image creation method capable of creating a more uniform image by suppressing a difference in spatial resolution and noise characteristics depending on a position inside a measurement site.
  • DESCRIPTION OF SYMBOLS 10 ... Biological measuring device, 14 ... Operation part, 16 ... Light emission / detection end, 22 ... Light source, 24 ... Optical switch, 26 ... Optical fiber for light source, 28 ... Optical fiber for emission, 30 ... Optical detector, 32 ... Shutter 34... Detection optical fiber 36... Signal processing circuit 38.

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Biomedical Technology (AREA)
  • Pathology (AREA)
  • Veterinary Medicine (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Biophysics (AREA)
  • Signal Processing (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radiology & Medical Imaging (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Psychiatry (AREA)
  • Physiology (AREA)
  • Oncology (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

　生体計測装置１０は、被計測部位Ｂに光を照射する光照射部と、被計測部位からの拡散光を検出する光検出部と、被計測部位の内部に関する再構成画像を作成する演算部１４とを備える。演算部１４は、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（Ｊは再構成画像の画素数）を算出し、次の反復式（但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより再構成画像を作成する。これにより、被計測部位内部の位置による空間分解能や雑音特性の差を抑制し、より均一に近い画像を作成することができる生体計測装置および画像作成方法が実現される。

Description

生体計測装置および画像作成方法
　本発明は、生体計測装置および画像作成方法に関するものである。
　頭部や乳房といった生体の内部情報を非侵襲的に計測する装置として、生体の光吸収特性を利用して内部情報を得る、いわゆる拡散光トモグラフィ（ＤＯＴ；Diffuse Optical
Tomography）を用いたものが提案されている。このような計測装置においては、計測対象となる生体の部位に対して所定の照射位置から光を照射し、当該部位の内部を散乱されつつ伝播された光を所定の検出位置で検出し、その強度や時間波形などの測定結果から、当該部位の内部情報、すなわち当該部位の内部に存在する腫瘍などの光吸収体に関する情報を得ることができる。なお、特許文献１には、拡散光トモグラフィを用いた生体計測方法が記載されている。また、非特許文献１及び２には、拡散光トモグラフィにおける逐次近似画像再構成法が記載されている。
特開２００１－２６４２４５号公報
Y. Ueda, K. Ohta, M. Oda, M. Miwa, Y. Tsuchiya, and Y. Yamashita,"Three-dimensional imaging of a tissuelike phantom by diffusion opticaltomography", Applied Optics Vol.40 No.34, pp.6349-6355 (2001) Y. Ueda, T. Yamanaka, D. Yamashita, T. Suzuki, E. Ohmae, M. Oda andY. Yamashita, "Reflectance Diffuse Optical Tomography: Its Application toHuman Brain Mapping", Japanese Journal of Applied Physics Vol.44 No.38, pp.L1203-L1206(2005)
　拡散光トモグラフィによる生体計測では、被計測部位内部における位置によって空間分解能や雑音特性が異なるため、不均一な画像が生成される。図１７は、そのような現象を説明するための模式図であり、被計測部位１００と、被計測部位１００の表面に設けられた光照射部１０１及び光検出部１０２が示されている。被計測部位１００の内部では、光照射位置１０１から出射されて光検出位置１０２に達する光子の飛行時間が短いほど、飛行距離が短く飛行経路が限定される。逆に、光子の飛行時間が長いほど、飛行距離が長くなり飛行経路は限定されない。そして、光子の飛行時間が短いデータには、図１７に示される領域Ａ１、すなわち被計測部位の表面に近い領域を通過する飛行経路Ｒ１が多く含まれる。また、光子の飛行時間が長いデータには、図１７に示される領域Ａ２、すなわち被計測部位の表面から遠い領域を通過する飛行経路Ｒ２が多く含まれる。したがって、被計測部位の表面から遠い領域の情報量は、被計測部位の表面に近い領域の情報量より少なくなってしまう。これにより、被計測部位の表面に遠い領域の空間分解能や雑音が、被計測部位の表面から近い領域の空間分解能や雑音よりも大きくなってしまう。
　本発明は、このような問題点に鑑みてなされたものであり、被計測部位内部の位置による空間分解能や雑音特性の差を抑制し、より均一に近い画像を作成することができる生体計測装置および画像作成方法を提供することを目的とする。
　上述した課題を解決するために、本発明に係る第１の生体計測装置は、被検者の被計測部位に光を照射する光照射部と、被計測部位からの拡散光を検出する光検出部と、光検出部からの出力信号に基づいて、被計測部位内の光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備え、演算部は、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、次の反復式
Figure JPOXMLDOC01-appb-M000007
 （但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。
　また、本発明に係る第２の生体計測装置は、被検者の被計測部位に放射線若しくは音波を照射する照射部と、被計測部位からの拡散した放射線若しくは音波を検出する検出部と、検出部からの出力信号に基づいて、被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備え、演算部は、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（１）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。なお、本発明において、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。
　これらの生体計測装置では、再構成画像の各画素毎に設定されるＪ個の係数ｗ1～ｗJを用いて、画像再構成のための逐次近似演算を行っている。例えば、これらの係数ｗ1～ｗJを、Ｎ回の反復演算において収束速度が最も遅い領域に当該画素の収束速度を一致させるように設定することで、収束速度が均一化され、被計測部位内部の位置による空間分解能や雑音特性の差を抑制してより均一に近い画像を作成することが可能となる。
　また、生体計測装置は、演算部が、再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の収束速度（以下、最低収束速度という）ＣNを求め、０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、１からＭまでの各ｍについて、次の反復式
Figure JPOXMLDOC01-appb-M000008
 を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)～ｘJ (N)を算出し、該画素値ｘ1 (N)～ｘJ (N)から得られる各部分領域の収束速度が最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の係数ｗjとする構成としてもよい。これにより、上述した効果をより好適に得ることができる。
　また、生体計測装置は、演算部が、Ｎ回の反復演算において収束速度が最も遅い部分領域の係数ｗjを１とする構成としてもよい。これにより、上述した効果をより好適に得ることができる。
　本発明に係る第１の画像作成方法は、被検者の被計測部位に光を照射し、被計測部位からの拡散光を検出し、該検出信号に基づいて被計測部位内の光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する方法であって、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、次の反復式
Figure JPOXMLDOC01-appb-M000009
 （但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。
　また、本発明に係る第２の画像作成方法は、被検者の被計測部位に放射線若しくは音波を照射し、被計測部位からの拡散した放射線若しくは音波を検出し、該検出信号に基づいて被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する方法であって、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（３）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。なお、本発明において、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。
　これらの画像作成方法では、再構成画像の各画素毎に設定されるＪ個の係数ｗ1～ｗJを用いて、画像再構成のための逐次近似演算を行っている。例えば、これらの係数ｗ1～ｗJを、Ｎ回の反復演算において収束速度が最も遅い領域に当該画素の収束速度を一致させるように設定することで、収束速度が均一化され、被計測部位内部の位置による空間分解能や雑音特性の差を抑制してより均一に近い画像を作成することが可能となる。
　また、画像作成方法は、再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の収束速度（以下、最低収束速度という）ＣNを求め、０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、１からＭまでの各ｍについて、次の反復式
Figure JPOXMLDOC01-appb-M000010
 を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)～ｘJ (N)を算出し、該画素値ｘ1 (N)～ｘJ (N)から得られる各部分領域の収束速度が最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の係数ｗjとする構成としてもよい。これにより、上述した効果をより好適に得ることができる。
　また、画像作成方法は、Ｎ回の反復演算において収束速度が最も遅い部分領域の係数ｗjを１とする構成としてもよい。これにより、上述した効果をより好適に得ることができる。
　本発明による生体計測装置および画像作成方法によれば、被計測部位内部の位置による空間分解能や雑音特性の差を抑制し、より均一に近い画像を作成することができる。
図１は、一実施形態に係る生体計測装置の構成を示す図である。 図２は、係数（ステップサイズ）の具体的な決定方法を示すフローチャートである。 図３は、シミュレーションにおいて再構成の対象となる２種類の被測定部位を示す図である。 図４は、被計測部位の内部における光子の飛行をシミュレーションする順問題解析における条件を示す図表である。 図５は、検出された光子ヒストグラムから画像を再構成する逆問題解析における条件を示す図表である。 図６は、シミュレーションにおいてステップサイズを決定するために使用された４つの画像を示す図である。 図７は、（ａ）図３（ａ）に示された画像に対応する、シミュレーションによる再構成後の画像、及び（ｂ）図３（ａ）に示された画像に対応する、ステップサイズを用いない従来の方法による再構成後の画像である。 図８は、図９に示される３本のライン上の画素値の変化を示すグラフである。 図９は、画像上に想定された３本のラインを示す図である。 図１０は、（ａ）図３（ｂ）に示された画像に対応する、シミュレーションによる再構成後の画像、及び（ｂ）図３（ｂ）に示された画像に対応する、ステップサイズを用いない従来の方法による再構成後の画像である。 図１１は、図１２に示される２本のライン上の画素値の変化を示すグラフである。 図１２は、画像上に想定された２本のラインを示す図である。 図１３は、（ａ）図３（ａ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像、及び（ｂ）図３（ａ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。 図１４は、（ａ）図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像、及び（ｂ）図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。 図１５は、部分領域毎のＮＳＤを示すグラフである。 図１６は、標準偏差の算出結果を示す図表である。 図１７は、拡散光トモグラフィにおいて不均一な画像が生成される現象を説明するための模式図である。
　以下、添付図面を参照しながら本発明による生体計測装置および画像作成方法の実施の形態を詳細に説明する。なお、図面の説明において同一の要素には同一の符号を付し、重複する説明を省略する。
　図１は、本発明の一実施形態に係る生体計測装置１０の構成を示す図である。本実施形態の生体計測装置１０は、計測対象である被検者の被計測部位Ｂに光を照射し、拡散光（戻り光）を検出し、その検出位置と計測された光量データ（例えば時間分解光子ヒストグラム）とに基づいて、光子の平均飛行経路と平均光路長を推定し、画像再構成問題として体内の情報を画像化する装置である。この装置によって得られる画像は、例えば腫瘍の位置や酸素化ヘモグロビン及び脱酸素化ヘモグロビンの分布を可視化したものであり、体組織の機能画像である。なお、被計測部位Ｂとしては、例えば頭部や女性の乳房等が想定される。
　生体計測装置１０は、計測光を被計測部位Ｂの内部に照射する光照射部と、光照射部からの光の照射により被計測部位Ｂから生じた拡散光を検出する光検出部と、光検出部からの出力信号に基づいて被計測部位Ｂの吸収係数の空間的分布を計算し、被計測部位Ｂの再構成画像を作成する演算部１４とを備えている。
　本実施形態の光照射部は、被計測部位Ｂに取り付けられるｎ個の光出射／検出端１６それぞれが有する光出射端、光源２２、および光スイッチ２４によって構成されている。光源２２としては、例えばレーザダイオードを使用することができる。計測光の波長としては、生体の透過率と定量すべき吸収体の分吸収係数との関係等から、７００ｎｍ～９００ｎｍ程度の近赤外線領域の波長が好ましい。
　計測光は、例えば連続光として光源２２から出射される。光源２２から出射された計測光は、光出射／検出端１６から被計測部位Ｂへ照射される。光スイッチ２４は、１入力ｎ出力の光スイッチであり、光源２２から光源用光ファイバ２６を介して光を入力し、この光を上記ｎ個の光出射／検出端１６それぞれに対して順に供給する。すなわち、光スイッチ２４は、各光出射／検出端１６に接続されたｎ本の出射用光ファイバ２８を１本ずつ順に選択し、当該出射用光ファイバ２８と光源２２とを光学的に接続する。
　本実施形態の光検出部は、前述したｎ個の光出射／検出端１６それぞれが有する光検出端と、ｎ個の光出射／検出端１６それぞれに対応するｎ個の光検出器３０と、各光検出器の入力部前段に配置されたｎ個のシャッター３２とによって構成されている。ｎ個の光検出器３０それぞれには、各光出射／検出端１６の光検出端に入射した被計測部位Ｂからの拡散光が、検出用光ファイバ３４を介して入力される。光検出器３０は、対応する光出射／検出端１６に到達した拡散光の光強度に応じてアナログ信号を生成する。光検出器３０としては、光電子増倍管（ＰＭＴ：Photomultiplier Tube）の他、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、様々なものを使用することができる。被計測部位Ｂからの拡散光が微弱であるときは、高感度あるいは高利得の光検出器を使用することが好ましい。光検出器３０の信号出力端には信号処理回路３６が接続されており、信号処理回路３６は、光検出器３０から出力されたアナログ信号をＡ／Ｄ変換して拡散光の光強度に応じたディジタル信号を生成し、該ディジタル信号を演算部１４へ提供する。
　演算部１４は、信号処理回路３６から提供されたディジタル信号に基づいて、被計測部位Ｂ内の光吸収係数分布を演算し、被計測部位Ｂの内部に関する再構成画像を作成する。演算部１４は、例えば、ＣＰＵ（Central Processing
Unit）といった演算手段及びメモリなどの記憶手段を有するコンピュータによって実現される。なお、演算部１４は、光源２２の発光、光スイッチ２４の動作及びシャッター３２の開閉を制御する機能を更に有するとよい。また、演算部１４には記録／表示部３８が接続されており、演算部１４における演算結果、すなわち被計測部位Ｂの再構成画像を可視化することが可能となっている。
　被計測部位Ｂの内部情報の算出すなわち内部情報計測は、例えば次のようにして行われる。ｎ個の光出射／検出端１６のそれぞれから被計測部位Ｂの内部へ計測光を順に照射し、被計測部位Ｂを通過して拡散した光を、ｎ個の光出射／検出端１６を介してｎ個の光検出器３０により検出する。この検出結果に基づいて、被計測部位Ｂの内部における吸収係数の空間的分布を演算し、腫瘍などの吸収体の位置や形状に関する情報（内部情報）を含む再構成画像を作成する。
　なお、演算部１４における吸収係数分布の算出には、例えば特許文献１に詳しく説明されているような周知の方法を用いるとよい。
　続いて、光吸収係数の空間的分布に基づいて、再構成画像を作成する方法について説明する。なお、以下に説明する演算は、演算部１４において行われる。ここで、拡散光トモグラフィにおける画像再構成問題を定式化するために、未知の光吸収係数分布に基づく再構成画像を構成する各画素の値を、次のＪ次元列ベクトルｘによって表す。
　　ｘ＝（ｘ1，ｘ2，…，ｘJ）T
また、光検出部において検出された測定データである光子ヒストグラムを、次のＩ次元列ベクトルＴによって表す。
　　Ｔ＝（Ｔ1，Ｔ2，…，ＴI）T
また、これらｘとＴとを相互に関連づけるＩ×Ｊ型のシステム行列Ｌを
　　Ｌ＝｛ｌij｝
と定義する。ｌiをＬのｉ行目の要素ベクトルとする。更に、光吸収係数分布が一様で既知の画像を、Ｊ次元列ベクトルｘrefとし、そのｘrefに対応する測定データである光子ヒストグラムを、次のＩ次元列ベクトルＢによって表す。
　　Ｂ＝（Ｂ1，Ｂ2，…，ＢI）T
　測定データである光子ヒストグラムに統計雑音が含まれていない場合、次式（５）が成り立つ。
Figure JPOXMLDOC01-appb-M000011
 しかし、統計雑音が混入した場合、上式（５）が成り立たない。したがって、統計雑音が混入している状態での最適なｘを求める必要がある。そこで、本実施形態ではいわゆる最尤推定法を用いてこのようなｘを求める。最尤推定法では、光検出部における光子の検出確率から尤度関数を立式し、その尤度関数を目的関数として最適化問題を解くことによって最適なｘを求めることができる。
　拡散光トモグラフィにおける光子の検出確率はポアソン分布に従い、その統計雑音もポアソン分布に従う。したがって、拡散光トモグラフィの最適化問題は、次式（６）によって表される。
Figure JPOXMLDOC01-appb-M000012
 また、式（６）における目的関数F（ｘ）は、次式（７）に示される対数尤度関数によって表される。
Figure JPOXMLDOC01-appb-M000013
　本実施形態では、ＯＳ－Ｃｏｎｖｅｘ法を用いて上式（６）の最大化問題を勾配法で解くことにより、画像再構成を行う。ＯＳ－Ｃｏｎｖｅｘ法では、測定データをサブセットと呼ばれる部分データ集合に分割し、各サブセットに対応する評価関数が増加するようにサブセット毎に解を更新する。これにより、１回の反復演算においてサブセットの数だけ解を更新できるので、収束速度が向上する。
　ｋを１からＮまでの整数（Ｎは反復演算の回数）とし、サブセット数をＱとし、ｑを１からＱまでの整数とし、第ｑ番目のサブセットのデータ集合をＳｑと定義する。具体的なＯＳ－Ｃｏｎｖｅｘ法の反復式は、次式（８）及び（９）によって表される。
Figure JPOXMLDOC01-appb-M000014
Figure JPOXMLDOC01-appb-M000015
 ここで、一般的な逐次近似解法の反復式は、次式（１０）に示される。この反復式によれば、第（ｋ＋１）番目の値ｘj (k+1)は、第ｋ番目の値ｘj (k)と更新量ｄj (k)との和によって表される。
Figure JPOXMLDOC01-appb-M000016
 なお、上述したＯＳ－Ｃｏｎｖｅｘ法では、更新量ｄj (k)は次式（１１）のように表される。
Figure JPOXMLDOC01-appb-M000017
　解の収束（収束速度）が早い画素では、Ｎ回の反復演算それぞれにおける更新量ｄj (1)～ｄj (N)の値が大きい。一方、収束速度が遅い画素では、Ｎ回の反復演算それぞれにおける更新量ｄj (1)～ｄj (N)の値が小さい。したがって、本来は同値となるべき２つの画素の収束速度が互いに異なると、一定回数の反復演算後におけるそれぞれの値が互いに異なる結果となる。つまり、画素毎の収束速度が大きく異なると、再構成画像における空間分解能の均一性が損なわれてしまうこととなる。
　すなわち、Ｊ個の画素における値ｘ1 (N)～ｘJ (N)の収束速度を互いに近づけることにより、Ｎ回の反復演算後における再構成画像の空間分解能（画質）を均一化することができる。その為には、値ｘ1 (N)～ｘJ (N)の収束速度が互いに近づく（好ましくは、略等しくなる）ように、上式（１０）における更新量ｄj (k)に対し画素毎に異なる任意の係数を乗ずるとよい。但し、更新量ｄj (k)は、当該第ｊ番目の画素における評価関数を増加させる最大の更新量であるので、収束速度が遅い画素の更新量ｄj (k)に対して１より大きい係数を乗ずると、評価関数は増加せず減少してしまうので好ましくない。そこで、収束速度が早い画素の更新量ｄj (k)に対し、１より小さい係数を乗算することによって、評価関数を増加させつつ、収束速度が遅い画素に合わせて収束速度を揃えることが可能となる。換言すれば、次式（１２）を満たす係数ｗj（以下、ステップサイズという）を各画素に設定し、このステップサイズｗjを更新量ｄj (k)に乗ずることによって、各画素値ｘ1 (N)～ｘJ (N)の収束速度を任意に制御することが可能となる。
Figure JPOXMLDOC01-appb-M000018
 したがって、上述した反復演算式（１０）は、次式（１３）のように書き換えられる。
Figure JPOXMLDOC01-appb-M000019
 例えば、上式（１３）に示される反復演算式を前述したＯＳ－Ｃｏｎｖｅｘ法に適用すると、式（９）に示された反復演算式は次式（１４）のようになる。
Figure JPOXMLDOC01-appb-M000020
　続いて、ステップサイズｗjを決定するための方法について説明する。
　＜収束速度の評価＞
　ステップサイズｗjを決定するためには、測定データに基づいて画素毎の収束速度を求めたのち、収束速度が最も遅い画素を判定する必要がある。画素毎の収束速度は、コントラスト回復率（ＣＲＣ：Contrast Recovery
Coefficient）によって評価することができる。
　拡散光トモグラフィの性格として、相互に隣接する画素の収束速度は互いに近いと考えられる。したがって本実施形態では、複数の画素をそれぞれ含む複数の部分領域に画像を分割し、スポットとなる領域を各部分領域内に配置し、該スポットにおけるＣＲＣをその領域の収束速度に相当する値として用いる。なお、ＣＲＣは、次式（１５）によって定義される。
Figure JPOXMLDOC01-appb-M000021
 但し、ＳＰはスポット領域における画素値を表し、ＢＧはスポット領域外（背景領域）の画素値を表す。また、添字ｍは領域内の平均値であることを示し、添字Ｒは画像再構成後の画素値であることを示し、添字Ｔｒは測定データに基づく画素値であることを示す。
　＜ステップサイズｗjの決定＞
　図２は、上述した収束速度の評価方法に基づく、ステップサイズｗjの具体的な決定方法を示すフローチャートである。この方法では、まず、画像をＥ個の部分領域に分割する（ステップＳ１）。なお、以下の説明において、或る任意の部分領域ｅの画素値の集合をＲｅとする。
　次に、任意の反復回数Ｎを設定し、Ｎ回の反復演算においてＣＲＣが最も小さくなる（すなわち、収束速度が最も遅い）部分領域を判定し、そのＣＲＣを最低収束速度ＣNと定義する（ステップＳ２）。そして、０＜ｖm＜１を満たすＭ個の任意の数値ｖm（但し、ｍは１からＭまでの整数であり、ｖ1～ｖMは互いに異なる値である）を用意する（ステップＳ３）。なお、このステップＳ３は、ステップＳ２やステップＳ１より前に行ってもよい。
　続いて、１からＭまでの各ｍについて、次の反復式（１６）を用いてＮ回の反復演算を行うことにより、Ｊ個の各画素の画素値ｘ1 (N)～ｘJ (N)を算出し、再構成画像を作成する（ステップＳ４）。
Figure JPOXMLDOC01-appb-M000022
 そして、該画素値ｘ1 (N)～ｘJ (N)から得られる或る部分領域の収束速度が最低収束速度ＣNと略一致した場合に、そのときの数値ｖmを該部分領域に含まれる複数の画素のステップサイズｗjとする（ステップＳ５）。以降、前述したステップＳ４及びＳ５をＭ回繰り返すか、或いは、収束速度が最も遅い部分領域を除く全ての部分領域についてステップサイズｗjが決定するまで、ステップＳ４及びＳ５を繰り返す。
　最後に、Ｎ回の反復演算において収束速度が最も遅い部分領域のステップサイズｗjを１とする（ステップＳ６）。こうして、全ての部分領域の全ての画素について、ステップサイズｗjが決定される。
　以上に説明した本実施形態に係る生体計測装置１０及び画像作成方法によって得られる効果について、従来方法の課題とともに説明する。
　拡散光トモグラフィは、生体内での近赤外光の拡散的な飛行経路を用いた画像再構成法である。拡散光トモグラフィの画像再構成では、解析的手法ではなく逐次近似画像再構成法が用いられるが、共役勾配法といった通常の逐次近似画像再構成法では、再構成画像内での空間分解能や雑音特性のばらつきが大きい不均一な画像が生成される問題がある。
　画質の不均一性の原因は、拡散光トモグラフィにおいては測定データに含まれる情報量が被測定部位の内部の位置によって大きく異なるためである。すなわち、時間分解計測法を用いた拡散光トモグラフィーでは、被計測部位表面の光入射端から入射された光子は被計測部位内で散乱を繰り返しながら飛行し、被計測部位表面の光検出端に到達して検出される。この入射から検出までの光子の飛行時間が短いほど、飛行距離が短くなり飛行経路も限られる。逆に、飛行時間が長いほど、飛行距離が長くなり飛行経路は限定されない。このため、測定データに含まれる情報量は、光子の飛行時間に応じて変化する。典型的には、光子の飛行時間が短い場合には、被計測部位表面に近い飛行経路の割合が多くなり、また、飛行時間が長い場合には、被計測部位表面から遠い（被計測部位の中央付近の）飛行経路の割合が多くなる。そのため、被計測部位の中央付近の情報量は表面付近と比べて少なくなってしまう。これにより、再構成画像の周辺部における空間分解能及び雑音が再構成画像の中央付近における空間分解能及び雑音と比べて大きくなるといった、再構成画像内での空間分解能及び雑音のばらつきが生じてしまう。
　本実施形態の生体計測装置１０及び画像作成方法では、再構成画像の各画素毎に設定されるＪ個の係数ｗ1～ｗJを用いて、画像再構成のための逐次近似演算を行う。そして、本実施形態のように、これらの係数ｗ1～ｗJを、Ｎ回の反復演算において収束速度が最も遅い部分領域に当該画素の収束速度を一致させるように設定することで、収束速度が均一化され、被計測部位内部の位置による空間分解能や雑音特性の差を抑制して、より均一に近い画像を作成することが可能となる。
　ここで、本実施形態に係る画像作成方法による上記効果を確認するために行われたシミュレーションの結果について説明する。図３（ａ）及び図３（ｂ）は、このシミュレーションにおいて再構成の対象となる２種類の被測定部位を示す図であり、各図には、例えば腫瘍といった複数の光吸収体Ｄが示されている。なお、図３（ａ）に示される３つの光吸収体Ｄの中心座標（ｘ，ｙ）は、それぞれ（２１，２０）、（６４，１０６）、及び（７５，５１）である。また、図３（ｂ）に示される４つの光吸収体Ｄの中心座標（ｘ，ｙ）は、それぞれ（５６，７９）、（６６，２４）、（１０９，４９）、及び（１０９，１１０）である。
　このシミュレーションでは、次のステップ（１）～（３）を行った。
（１）各画素の収束速度を評価する為に、画像内において均等に配置されたスポットにおける画像再構成を行い、ステップサイズｗ1～ｗJを決定した。
（２）得られたステップサイズｗ1～ｗJを用いて画像の再構成を行い、再構成画像の全体で空間分解能が略均一となっていることを確認した。
（３）背景画像の再構成による領域間での雑音伝播を評価した。
　なお、上記ステップ（１）～（３）に共通する条件は、次のとおりである。まず、被計測部位の内部における光子の飛行をシミュレーションする順問題解析（Forward problem）と、検出された光子ヒストグラムから画像を再構成する逆問題解析（Inverse problem）という２種類の解析における条件を、それぞれ図４及び図５に示す。図４及び図５に示されるように、順問題解析と逆問題解析との間でグリッドサイズや画像サイズが異なるのは、連続系を離散系に変換したことによる。光検出部において検出される光子は、必ず離散的なデータとして取得される。一方、数値シミュレーションの値は連続関数に近い形で計算される。したがって、本シミュレーションでは、ダウンサンプリングすることにより連続系を離散系に変換している。
　（１）ステップサイズｗ1～ｗJの決定
　測定データに基づく１３２行１３２列の画素を含む画像を１２行１２列の部分領域に分割し、各部分領域の中央にホットスポットを配置した。本シミュレーションでは、図６（ａ）～図６（ｄ）に示される４つの画像を再構成することにより、各部分領域のＣＲＣを求め、前述した方法によって全画素のステップサイズｗ1～ｗJを決定した。
　（２）画像の再構成
　次に、ステップ（１）により決定されたステップサイズｗ1～ｗJを用い、式（１４）に示された反復演算式を演算することにより、画像の再構成を行った。なお、このステップでは、各部分領域のＣＲＣに差があると画像中にエッジとなって現れてしまうので、ステップサイズｗ1～ｗJを各画素値とする画像（ステップサイズ画像）を平均値フィルターを用いてぼかし、滑らかに変化するステップサイズ画像を得た。本シミュレーションでは、カーネルサイズ９×９の平均値フィルターを用いた。
　図７（ａ）は、図３（ａ）に示された画像に対応する、本シミュレーションによる再構成後の画像である。図７（ｂ）は、図３（ａ）に示された画像に対応する、ステップサイズｗ1～ｗJを用いない従来の方法による再構成後の画像である。図８（ａ）～図８（ｃ）は、図９に示されるように光吸収体Ｄを通る３本のラインＬ１～Ｌ３上の画素値の変化を示すグラフであって、グラフＧ２１は本シミュレーションによる画素値の変化を示し、グラフＧ２２は従来の方法による画素値の変化を示している。また、図１０（ａ）は、図３（ｂ）に示された画像に対応する、本シミュレーションによる再構成後の画像である。図１０（ｂ）は、図３（ｂ）に示された画像に対応する、ステップサイズｗ1～ｗJを用いない従来の方法による再構成後の画像である。図１１（ａ）及び図１１（ｂ）は、図１２に示されるように光吸収体Ｄを通る２本のラインＬ４及びＬ５上の画素値の変化を示すグラフであって、グラフＧ２１は本シミュレーションによる画素値の変化を示し、グラフＧ２２は従来の方法による画素値の変化を示している。
　図７（ｂ）を参照すると、従来の再構成画像では、左上といった画像周辺部の光吸収体Ｄの画素値が比較的高くなっており、画像中央付近の光吸収体Ｄの画素値と大きく異なっていることがわかる。これに対し、ステップサイズｗ1～ｗJを用いた本実施形態による再構成画像では、図７（ａ）に示されるように、各光吸収体Ｄの画素値が、画像の周辺部及び中央付近においてほぼ均一な値となっている。また、図８を参照すると、従来の再構成画像では、画像周辺部における画素値の変化が大きい。これに対し、図８（ａ）～図８（ｃ）に示されるように、本実施形態による再構成画像では画像周辺部における画素値の変化が小さくなっている。これは、本実施形態において画像中央付近ではステップサイズｗjが大きくなって反復演算における更新量が大きくなり、画像周辺部ではステップサイズｗjが小さくなって更新量が小さくなったことによる。以上の結果から、本実施形態の生体計測装置１０及び画像作成方法によれば、収束速度が均一化されたことにより、被計測部位内部の位置による空間分解能の差を抑制して、より均一に近い画像を作成できることがわかる。
　次に、統計雑音を付加した場合のシミュレーション結果を示す。このシミュレーションでは、各検出光子ヒストグラムの最大値が一定値（例えば５０）となるように光子ヒストグラムをフィッティングした後、この光子ヒストグラムにポアソン雑音を付加した。図１３（ａ）は、図３（ａ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。図１３（ｂ）は、図３（ａ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。また、図１４（ａ）は、図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。図１４（ｂ）は、図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。
　図１３（ｂ）及び図１４（ｂ）を参照すると、統計雑音が存在する場合、ステップサイズｗ1～ｗJを用いない従来の方法では、画像周辺部の空間分解能が高くなって雑音の影響が大きく出ているが、画像中央付近では空間分解能が高くなく雑音の影響が小さいことがわかる。これに対し、図１３（ａ）及び図１４（ａ）を参照すると、ステップサイズｗ1～ｗJを用いた本実施形態による方法では、反復演算における画像周辺部の更新量が抑えられているので、雑音の影響が顕著に抑制され、左上に存在する光吸収体Ｄの形を明確に視認できる。一方、画像中央付近では、雑音の影響が従来方法と同程度に抑えられている。なお、画像中央付近の統計雑音が画像周辺部の統計雑音よりも大きいのは、統計雑音の有無によって分解能の向上速度に変化が生じるためであり、統計雑音がない場合を仮定して決定したステップサイズを用いたことに起因する現象である。
　このように、計測データに統計雑音が含まれる場合、本実施形態の生体計測方法では従来方法と比較して画素値に変化はあるものの、画像全体の統計雑音の形態的特徴に大きな変化がないことがわかる。すなわち、本実施形態の生体計測方法によれば、画像周辺部において選択的に、分解能の向上を遅くすることができていることがわかる。
　続いて、光吸収体Ｄを配置しない背景画像に関する雑音評価について説明する。ここでは、光吸収体Ｄが存在しない計測データにポアソン雑音を付加した。このときの再構成画像は、雑音によって歪んだ画像となる。この再構成画像を６行６列の３６個の部分領域に分割し、左上から右下へ順に番号を付与した。その後、それぞれの部分領域を雑音指標（ＮＳＤ：Normalized Standard
Deviation）によって評価した。なお、ＮＳＤの計算式は次式（１７）である。また、図１５は、部分領域毎のＮＳＤを示すグラフであって、横軸は部分領域の番号を示し、縦軸はＮＳＤの値を示している。なお、添字ＳＤは標準偏差であることを示す。
Figure JPOXMLDOC01-appb-M000023
 また、部分領域毎のＮＳＤを母集団とする標準偏差を算出し、この標準偏差を部分領域間の雑音の不均一性を評価する値として用いる。図１６は、標準偏差の算出結果を示す図表である。
　図１５を参照すると、従来方法では部分領域間のＮＳＤのばらつきが大きい。これに対し、本実施形態の方法では、ＮＳＤのばらつきが全体的に抑えられていることがわかる。このことから、本実施形態の方法では、空間分解能を均一に近づけるために用いられるステップサイズｗjが、雑音伝播の不均一をも緩和していることがわかる。更に、図１６に示されるように、雑音の不均一性についても、本実施形態の方法では従来方法よりも小さくなる。
　本発明による生体計測装置および画像作成方法は、上述した実施形態に限られるものではなく、他に様々な変形が可能である。例えば、上記実施形態では被検者の被計測部位に光を照射し、被計測部位からの拡散光を検出し、該検出信号に基づいて被計測部位内における光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成している。しかしながら、本発明は、光を用いた生体計測装置や画像作成方法に限らず、放射線や音波を用いた生体計測装置や画像作成方法にも適用可能である。
　すなわち、本発明に係る生体計測装置は、被検者の被計測部位に放射線若しくは音波を照射する照射部と、被計測部位からの拡散した放射線若しくは音波を検出する検出部と、検出部からの出力信号に基づいて、被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備えてもよい。この場合、演算部は、上述した実施形態と同様の手法により、再構成画像を作成するとよい。また、本発明に係る画像作成方法は、被検者の被計測部位に放射線若しくは音波を照射し、被計測部位からの拡散した放射線若しくは音波を検出し、該検出信号に基づいて被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を上述した実施形態と同様の手法により作成する方法であってもよい。ここで、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。
　上記実施形態による第１の生体計測装置では、被検者の被計測部位に光を照射する光照射部と、被計測部位からの拡散光を検出する光検出部と、光検出部からの出力信号に基づいて、被計測部位内の光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備え、演算部は、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（１）を用いた逐次近似演算を行うことにより再構成画像を作成する構成としている。
　また、上記実施形態による第２の生体計測装置では、被検者の被計測部位に放射線若しくは音波を照射する照射部と、被計測部位からの拡散した放射線若しくは音波を検出する検出部と、検出部からの出力信号に基づいて、被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備え、演算部は、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（１）を用いた逐次近似演算を行うことにより再構成画像を作成する構成としている。なお、本構成において、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。
　また、生体計測装置は、演算部が、再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の収束速度（以下、最低収束速度という）ＣNを求め、０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、１からＭまでの各ｍについて、上記反復式（２）を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)～ｘJ (N)を算出し、該画素値ｘ1 (N)～ｘJ (N)から得られる各部分領域の収束速度が最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の係数ｗjとする構成としてもよい。これにより、上述した効果をより好適に得ることができる。
　また、生体計測装置は、演算部が、Ｎ回の反復演算において収束速度が最も遅い部分領域の係数ｗjを１とする構成としてもよい。これにより、上述した効果をより好適に得ることができる。
　上記実施形態による第１の画像作成方法では、被検者の被計測部位に光を照射し、被計測部位からの拡散光を検出し、該検出信号に基づいて被計測部位内の光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する方法であって、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（３）を用いた逐次近似演算を行うことにより再構成画像を作成する構成としている。
　また、上記実施形態による第２の画像作成方法では、被検者の被計測部位に放射線若しくは音波を照射し、被計測部位からの拡散した放射線若しくは音波を検出し、該検出信号に基づいて被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する方法であって、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（３）を用いた逐次近似演算を行うことにより再構成画像を作成する構成としている。なお、本構成において、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。
　また、画像作成方法は、再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の収束速度（以下、最低収束速度という）ＣNを求め、０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、１からＭまでの各ｍについて、上記反復式（４）を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)～ｘJ (N)を算出し、該画素値ｘ1 (N)～ｘJ (N)から得られる各部分領域の収束速度が最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の係数ｗjとする構成としてもよい。これにより、上述した効果をより好適に得ることができる。
　また、画像作成方法は、Ｎ回の反復演算において収束速度が最も遅い部分領域の係数ｗjを１とする構成としてもよい。これにより、上述した効果をより好適に得ることができる。
　本発明は、被計測部位内部の位置による空間分解能や雑音特性の差を抑制し、より均一に近い画像を作成することができる生体計測装置および画像作成方法として利用可能である。
　１０…生体計測装置、１４…演算部、１６…光出射／検出端、２２…光源、２４…光スイッチ、２６…光源用光ファイバ、２８…出射用光ファイバ、３０…光検出器、３２…シャッター、３４…検出用光ファイバ、３６…信号処理回路、３８…表示部。

Claims (8)

  1. 　被検者の被計測部位に光を照射する光照射部と、
    　前記被計測部位からの拡散光を検出する光検出部と、
    　前記光検出部からの出力信号に基づいて、前記被計測部位内の光吸収係数分布を演算し、前記被計測部位の内部に関する再構成画像を作成する演算部と
    を備え、
    　前記演算部は、前記再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは前記再構成画像の画素数である。）を算出し、次の反復式
    Figure JPOXMLDOC01-appb-M000001
     （但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、生体計測装置。
  2. 　被検者の被計測部位に放射線若しくは音波を照射する照射部と、
    　前記被計測部位からの拡散した前記放射線若しくは前記音波を検出する検出部と、
    　前記検出部からの出力信号に基づいて、前記被計測部位内における前記放射線若しくは前記音波の吸収係数分布を演算し、前記被計測部位の内部に関する再構成画像を作成する演算部と
    を備え、
    　前記演算部は、前記再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは前記再構成画像の画素数である。）を算出し、次の反復式
    Figure JPOXMLDOC01-appb-M000002
     （但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、生体計測装置。
  3. 　前記演算部は、
    　前記再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記収束速度（以下、最低収束速度という）ＣNを求め、
    　０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、
    　１からＭまでの各ｍについて、次の反復式
    Figure JPOXMLDOC01-appb-M000003
     を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)～ｘJ (N)を算出し、
    　該画素値ｘ1 (N)～ｘJ (N)から得られる各部分領域の前記収束速度が前記最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の前記係数ｗjとすることを特徴とする、請求項１または２に記載の生体計測装置。
  4. 　前記演算部は、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記係数ｗjを１とすることを特徴とする、請求項３に記載の生体計測装置。
  5. 　被検者の被計測部位に光を照射し、前記被計測部位からの拡散光を検出し、該検出信号に基づいて前記被計測部位内の光吸収係数分布を演算し、前記被計測部位の内部に関する再構成画像を作成する方法であって、
    　前記再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは前記再構成画像の画素数である。）を算出し、
    　次の反復式
    Figure JPOXMLDOC01-appb-M000004
     （但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、画像作成方法。
  6. 　被検者の被計測部位に放射線若しくは音波を照射し、前記被計測部位からの拡散した前記放射線若しくは前記音波を検出し、該検出信号に基づいて前記被計測部位内における前記放射線若しくは前記音波の吸収係数分布を演算し、前記被計測部位の内部に関する再構成画像を作成する方法であって、
    　前記再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは前記再構成画像の画素数である。）を算出し、
    　次の反復式
    Figure JPOXMLDOC01-appb-M000005
     （但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、画像作成方法。
  7. 　前記再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記収束速度（以下、最低収束速度という）ＣNを求め、
    　０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、
    　１からＭまでの各ｍについて、次の反復式
    Figure JPOXMLDOC01-appb-M000006
     を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)～ｘJ (N)を算出し、
    　該画素値ｘ1 (N)～ｘJ (N)から得られる各部分領域の前記収束速度が前記最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の前記係数ｗjとすることを特徴とする、請求項５または６に記載の画像作成方法。
  8. 　Ｎ回の反復演算において収束速度が最も遅い部分領域の前記係数ｗjを１とすることを特徴とする、請求項７に記載の画像作成方法。
PCT/JP2012/066552 2011-07-07 2012-06-28 生体計測装置および画像作成方法 WO2013005635A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
EP12807701.3A EP2730231B1 (en) 2011-07-07 2012-06-28 Biometric apparatus and image-generating method
US14/130,690 US9282932B2 (en) 2011-07-07 2012-06-28 Biometric apparatus and image-generating method
CN201280033086.7A CN103635148B (zh) 2011-07-07 2012-06-28 生物体测量装置以及图像制作方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2011-151086 2011-07-07
JP2011151086A JP5782314B2 (ja) 2011-07-07 2011-07-07 生体計測装置および画像作成方法

Publications (1)

Publication Number Publication Date
WO2013005635A1 true WO2013005635A1 (ja) 2013-01-10

Family

ID=47436990

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2012/066552 WO2013005635A1 (ja) 2011-07-07 2012-06-28 生体計測装置および画像作成方法

Country Status (5)

Country Link
US (1) US9282932B2 (ja)
EP (1) EP2730231B1 (ja)
JP (1) JP5782314B2 (ja)
CN (1) CN103635148B (ja)
WO (1) WO2013005635A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108398245A (zh) * 2018-05-28 2018-08-14 南昌大学 一种同心发散光源的检测方法与装置

Families Citing this family (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6595211B2 (ja) * 2015-05-11 2019-10-23 キヤノンメディカルシステムズ株式会社 核医学診断装置及び核医学画像処理装置
CN104997533B (zh) * 2015-06-23 2017-11-14 武汉超信电子工程有限公司 超声探头几何参数自动校正方法和装置
US9730649B1 (en) 2016-09-13 2017-08-15 Open Water Internet Inc. Optical imaging of diffuse medium
WO2018189815A1 (ja) 2017-04-11 2018-10-18 オリンパス株式会社 光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置
US10778911B2 (en) 2018-03-31 2020-09-15 Open Water Internet Inc. Optical transformation device for imaging
US10506181B2 (en) 2018-03-31 2019-12-10 Open Water Internet Inc. Device for optical imaging
US10778912B2 (en) 2018-03-31 2020-09-15 Open Water Internet Inc. System and device for optical transformation
US10966612B2 (en) 2018-06-14 2021-04-06 Open Water Internet Inc. Expanding beam optical element
US10962929B2 (en) 2018-09-14 2021-03-30 Open Water Internet Inc. Interference optics for optical imaging device
US10874370B2 (en) 2019-01-28 2020-12-29 Open Water Internet Inc. Pulse measurement in optical imaging
US10955406B2 (en) 2019-02-05 2021-03-23 Open Water Internet Inc. Diffuse optical imaging with multiple beams
US11320370B2 (en) 2019-06-26 2022-05-03 Open Water Internet Inc. Apparatus for directing optical and acoustic signals
US11581696B2 (en) 2019-08-14 2023-02-14 Open Water Internet Inc. Multi-channel laser
US11622686B2 (en) 2019-11-22 2023-04-11 Open Water Internet, Inc. Optical imaging with unshifted reference beam
CN111332287B (zh) * 2020-03-11 2021-01-22 合肥鑫晟光电科技有限公司 一种驾驶辅助方法及系统
US20210330266A1 (en) * 2020-04-24 2021-10-28 Hi Llc Systems and Methods for Noise Removal in an Optical Measurement System
US11819318B2 (en) 2020-04-27 2023-11-21 Open Water Internet Inc. Optical imaging from light coherence
US11259706B2 (en) 2020-05-19 2022-03-01 Open Water Internet Inc. Dual wavelength imaging and out of sample optical imaging
US11559208B2 (en) 2020-05-19 2023-01-24 Open Water Internet Inc. Imaging with scattering layer

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001264245A (ja) 2000-03-21 2001-09-26 Hamamatsu Photonics Kk 散乱吸収体内部の光路分布計算方法
JP2009516561A (ja) * 2005-11-23 2009-04-23 コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ 混濁媒質の光学的撮像の方法および装置
JP2010094500A (ja) * 2008-09-19 2010-04-30 Canon Inc 測定装置及び測定方法
JP2010514478A (ja) * 2006-12-22 2010-05-06 ワシントン・ユニバーシティ 拡散光トモグラフィのための高性能イメージングシステム及び関連した使用方法
WO2010143421A1 (ja) * 2009-06-10 2010-12-16 パナソニック株式会社 光融合型イメージング方法、光融合型イメージング装置、プログラムおよび集積回路
JP2011505954A (ja) * 2007-12-17 2011-03-03 コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ 混濁媒体の内部における不均質性の存在を検出するための方法及び、混濁媒体の内部を画像化するための装置

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4662222A (en) 1984-12-21 1987-05-05 Johnson Steven A Apparatus and method for acoustic imaging using inverse scattering techniques
US5694938A (en) * 1995-06-07 1997-12-09 The Regents Of The University Of California Methodology and apparatus for diffuse photon mimaging
US5905261A (en) 1995-12-01 1999-05-18 Schotland; John Carl Imaging system and method using direct reconstruction of scattered radiation
US8831709B2 (en) * 2004-09-24 2014-09-09 Softscan Healthcare Group Ltd. Method for 3-dimensional fluorescence tomographic imaging
EP2239560B1 (en) 2007-12-27 2020-04-22 Omron Corporation X-ray examining apparatus and x-ray examining method
TWI394490B (zh) 2008-09-10 2013-04-21 Omron Tateisi Electronics Co Ｘ射線檢查裝置及ｘ射線檢查方法
CN101509871B (zh) 2009-03-25 2011-09-28 华中科技大学 适于小动物的荧光分子层析成像方法
JP5600946B2 (ja) 2010-01-28 2014-10-08 株式会社島津製作所 断層撮影装置

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001264245A (ja) 2000-03-21 2001-09-26 Hamamatsu Photonics Kk 散乱吸収体内部の光路分布計算方法
JP2009516561A (ja) * 2005-11-23 2009-04-23 コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ 混濁媒質の光学的撮像の方法および装置
JP2010514478A (ja) * 2006-12-22 2010-05-06 ワシントン・ユニバーシティ 拡散光トモグラフィのための高性能イメージングシステム及び関連した使用方法
JP2011505954A (ja) * 2007-12-17 2011-03-03 コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ 混濁媒体の内部における不均質性の存在を検出するための方法及び、混濁媒体の内部を画像化するための装置
JP2010094500A (ja) * 2008-09-19 2010-04-30 Canon Inc 測定装置及び測定方法
WO2010143421A1 (ja) * 2009-06-10 2010-12-16 パナソニック株式会社 光融合型イメージング方法、光融合型イメージング装置、プログラムおよび集積回路

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Y. UEDA; K. OHTA; M. ODA; M. MIWA; Y. TSUCHIYA; Y. YAMASHITA: "Three-dimensional imaging of a tissuelike phantom by diffusion optical tomography", APPLIED OPTICS, vol. 40, no. 34, 2001, pages 6349 - 6355
Y. UEDA; T. YAMANAKA; D. YAMASHITA; T. SUZUKI; E. OHMAE; M. ODA; Y. YAMASHITA: "Reflectance Diffuse Optical Tomography: Its Application to Human Brain Mapping", JAPANESE JOURNAL OF APPLIED PHYSICS, vol. 44, no. 38, 2005, pages L1203 - L1206

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108398245A (zh) * 2018-05-28 2018-08-14 南昌大学 一种同心发散光源的检测方法与装置

Also Published As

Publication number Publication date
CN103635148A (zh) 2014-03-12
CN103635148B (zh) 2015-11-25
EP2730231A1 (en) 2014-05-14
EP2730231A4 (en) 2015-03-25
US20140135620A1 (en) 2014-05-15
JP5782314B2 (ja) 2015-09-24
US9282932B2 (en) 2016-03-15
EP2730231B1 (en) 2017-11-29
JP2013019696A (ja) 2013-01-31

Similar Documents

Publication Publication Date Title
JP5782314B2 (ja) 生体計測装置および画像作成方法
Lyons et al. Computational time-of-flight diffuse optical tomography
JP5975437B2 (ja) 生体計測装置の計測データ選択方法、生体計測装置の光出射位置決定方法、および生体計測装置
US8817257B2 (en) Method for reconstructing the optical properties of a medium using a combination of a plurality of mellin-laplace transforms of a magnitude comprising a time distribution of a received signal, and associated reconstruction system
JP5570674B2 (ja) 物体観測装置、物体観測方法、および記録媒体
US10088571B2 (en) Underwater sensing system
WO2012153769A1 (ja) 光トモグラフィ装置
US20070244395A1 (en) Systems and methods for multi-spectral bioluminescence tomography
EP2302362B1 (fr) Dispositif et procédé de reconstruction spatiale d'une cartographie de fluorescence
JP7424289B2 (ja) 情報処理装置、情報処理方法、情報処理システム、およびプログラム
CN103815929B (zh) 被检体信息获取装置
Zuo et al. Spectral crosstalk in photoacoustic computed tomography
JP5924658B2 (ja) 濃度定量装置、光吸収係数算出方法、等価散乱係数算出方法、濃度定量方法、光吸収係数の算出を行うプログラム及び濃度の算出を行うプログラム
US5746211A (en) Absorption imaging system and method using direct reconstruction of scattered radiation
US20130158926A1 (en) Method for reconstructing the optical properties of a medium with computing of a signal corrected as a function of a first modeling function for a reference medium and of a second distribution for a medium to be characterized, and associated reconstruction system
JP5420163B2 (ja) 生体計測装置
JP2008527321A (ja) 混濁媒質の光学特性を決定するための方法
Kalyanov et al. Time-multiplexing approach for fast time-domain near-infrared optical tomography combined with neural-network-enhanced image reconstruction
Okawa et al. Reduction of Poisson noise in measured time-resolved data for time-domain diffuse optical tomography
Soonthornsaratoon Gradient-based methods for quantitative photoacoustic tomography
Żołek et al. Analysis of estimation of optical properties of sub superficial structures in multi layered tissue model using distribution function method
Munoz-Barrutia et al. Sparse algebraic reconstruction for fluorescence mediated tomography
Peter Comparative study on 3D modelling of breast cancer using Nir-Fdot in COMSOL
Lin et al. Anatomical Modeling and Optimization of Speckle Contrast Optical Tomography
WO2018189815A1 (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: 12807701

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

REEP Request for entry into the european phase

Ref document number: 2012807701

Country of ref document: EP

WWE Wipo information: entry into national phase

Ref document number: 2012807701

Country of ref document: EP

WWE Wipo information: entry into national phase

Ref document number: 14130690

Country of ref document: US