WO2021095852A1 - 信号処理装置、信号処理方法および信号処理プログラム - Google Patents

信号処理装置、信号処理方法および信号処理プログラム Download PDF

Info

Publication number
WO2021095852A1
WO2021095852A1 PCT/JP2020/042459 JP2020042459W WO2021095852A1 WO 2021095852 A1 WO2021095852 A1 WO 2021095852A1 JP 2020042459 W JP2020042459 W JP 2020042459W WO 2021095852 A1 WO2021095852 A1 WO 2021095852A1
Authority
WO
WIPO (PCT)
Prior art keywords
phase error
point
sample
path
bulk phase
Prior art date
Application number
PCT/JP2020/042459
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 US17/776,178 priority Critical patent/US20220390368A1/en
Priority to JP2021556177A priority patent/JPWO2021095852A1/ja
Publication of WO2021095852A1 publication Critical patent/WO2021095852A1/ja

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/41—Refractivity; Phase-affecting properties, e.g. optical path length
    • G01N21/45—Refractivity; Phase-affecting properties, e.g. optical path length using interferometric methods; using Schlieren methods
    • G01N21/455—Schlieren methods, e.g. for gradient index determination; Shadowgraph
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02055—Reduction or prevention of errors; Testing; Calibration
    • G01B9/02062—Active error reduction, i.e. varying with time
    • G01B9/02064—Active error reduction, i.e. varying with time by particular adjustment of coherence gate, i.e. adjusting position of zero path difference in low coherence interferometry
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/0209—Low-coherence interferometers
    • G01B9/02091—Tomographic interferometers, e.g. based on optical coherence

Definitions

  • the present invention relates to signal processing devices, signal processing methods and signal processing programs.
  • the present application claims priority with respect to Japanese Patent Application No. 2019-206506 filed in Japan on November 14, 2019, the contents of which are incorporated herein by reference.
  • OCT optical coherence tomography
  • FD Fourier domain OCT
  • Non-Patent Document 1 describes a method of performing digital refocusing by performing phase correction on a three-dimensional OCT image (three-dimensional OCT volume) by using a phase difference method.
  • the phase difference method is a method of estimating a bulk phase error by calculating a differential phase field in two directions ((x, y) directions) in front of a three-dimensional OCT volume and accumulating the differential phase field. ..
  • Non-Patent Document 2 an attempt was made to reduce the error by performing a smoothing process on the estimated bulk phase error. Since this method lacks a theoretical basis for error reduction, it may rather be a factor of lowering the estimation accuracy or a factor of failure in estimating the bulk phase error.
  • the present invention has been made in view of the above points, and an object of the present invention is to provide a signal processing device, a signal processing method, and a signal processing program capable of improving the estimation accuracy of bulk phase error. ..
  • the present invention has been made to solve the above problems, and one aspect of the present invention is to set a phase gradient in a plane intersecting the light irradiation direction of an optical interference fault signal representing the state of a sample.
  • the phase gradient calculation unit that calculates for each sample point arranged on the plane and the multiple paths from the starting point, which is the sample point for which the bulk phase error has been determined, to the end point, which is the sample point for which the bulk phase error has not been determined.
  • the phase gradient for each sample point is integrated along each path to calculate the path-specific phase error at the end point, and the path-specific phase error is combined between the plurality of paths to form a bulk phase error at the end point.
  • It is a signal processing apparatus including a bulk phase error calculation unit for determining.
  • Another aspect of the present invention is the signal processing device of (1), in which the bulk phase error calculation unit calculates an average value between paths of the phase-specific phase error as the bulk phase error. May be good.
  • Another aspect of the present invention is the signal processing apparatus of (2), in which the bulk phase error calculation unit uses a weighting coefficient based on the path length of each of the plurality of paths as the average value. The weighted average value of the phase error for each path may be calculated.
  • Another aspect of the present invention is the signal processing device according to any one of (1) to (3), and each of the plurality of routes may be included within a predetermined range from the end point. ..
  • Another aspect of the present invention is the signal processing device according to any one of (1) to (4), wherein each of the plurality of paths has n (n is 1 or more and N-1 or less).
  • the sample points may be within the range equal to or less than the spatial resolution of the optical system.
  • Another aspect of the present invention is the signal processing device according to any one of (1) to (5), in which a path pattern indicating the plurality of paths is associated with the positional relationship between the start point and the end point.
  • a plurality of routes may be set in advance, and a plurality of routes indicated by a route pattern satisfying the positional relationship may be selected with the target sample point for which the bulk phase error is calculated as the end point.
  • Another aspect of the present invention is the signal processing apparatus of (6), in which one sample point on the plane is set as a reference point in the bulk phase error calculation unit, and a sample close to the reference point is used.
  • the order in which the bulk phase error is calculated is predetermined so that the points are prioritized, and when the sample points whose order is the mth time (m is an integer of 1 or more) are set as the target sample points, the order is described. May select the plurality of routes represented by a route pattern satisfying the positional relationship, starting from the sample point of the m-1th time.
  • Another aspect of the present invention is the signal processing device according to any one of (1) to (7), wherein the bulk phase error at the end point is projected on the end point and the plane in common.
  • a signal correction unit that subtracts from the phase of the signal value of the optical interference fault signal at the sample points arranged in the three-dimensional space may be provided.
  • Another aspect of the present invention is a method in a signal processing apparatus, in which a sample point in which a phase gradient in a plane intersecting the light irradiation direction of an optical interference fault signal representing a sample state is arranged on the plane is arranged.
  • the phase for each sample point is calculated for each of the phase gradient calculation process and the plurality of paths from the starting point where the bulk phase error is determined to the end point where the bulk phase error is undetermined.
  • Bulk phase error calculation process in which the gradient is integrated along each path to calculate the path-specific phase error at the end point, and the path-specific phase error is combined between the plurality of paths to determine the bulk phase error at the end point. It is a signal processing method having.
  • Another aspect of the present invention is a phase gradient calculation procedure for calculating the phase gradient in a plane intersecting the light irradiation direction of the optical interference fault signal representing the state of the sample for each sample point arranged on the plane. , Integrate the phase gradient for each sample point along each of the multiple paths from the starting point, which is the sample point where the bulk phase error has been determined, to the ending point, which is the sample point where the bulk phase error has not been determined.
  • a signal processing program for calculating the path-specific phase error at the end point, synthesizing the path-specific phase error between the plurality of paths, and executing the bulk phase error procedure for determining the bulk phase error at the end point. Is.
  • the estimation accuracy of the bulk phase error can be improved.
  • FIG. 1 is a configuration diagram showing an example of an optical interference tomographic meter 1 according to the present embodiment.
  • the optical coherence tomography 1 constitutes an observation system for observing the state of a sample using OCT.
  • the optical interference tomometer 1 irradiates the sample Sm with light, and acquires and acquires the interference light generated by interfering the reflected light reflected from the sample Sm with the reference light reflected by the reference mirror 40 (described later). It is a device that generates an image showing the surface of the sample Sm and the state inside the sample Sm from the interference light.
  • the object to be observed as the sample Sm may be, for example, a human or animal living body or a non-living body.
  • the living body can be the fundus, blood vessels, teeth, subcutaneous tissue, and the like.
  • the non-living organism may be an artificial structure such as an electronic part or a mechanical part, a natural structure such as a stone or a mineral, or a substance having no specific shape.
  • the optical interference tomometer 1 includes a light source 10, a beam splitter 20, collimators 30a, 30b, 50a, 50b, a reference mirror 40, a galvanometer mirror 60a, 60b, a spectroscope 70, and a signal processing device 100.
  • the beam splitter 20, the collimators 30a, 30b, 50a, 50b, the reference mirror 40, the galvanometer mirrors 60a, 60b, and the spectroscope 70 constitute an optical system called an interferometer.
  • the interferometer illustrated in FIG. 1 is a Michelson interferometer including an optical fiber F.
  • the light source 10, the spectroscope 70, the collimator 30a, and the collimator 50a are each connected to the beam splitter by using an optical fiber F.
  • the optical fiber F has a transmission band including a wavelength band of light emitted from the light source 10.
  • the optical coherence tomography 1 is, for example, a Fourier domain OCT (FD-OCT: Fourier-Domain OCT).
  • the FD-OCT may be, for example, any of Spectral Domain OCT (SD-OCT), Wavelength Sweep OCT (Swept-Source OCT), and the like, but is not limited thereto.
  • the optical coherence tomography 1 may be, for example, a time domain OCT (Time-Domain OCT: TD-OCT).
  • TD-OCT Time-Domain OCT
  • the light source 10 irradiates, for example, probe light having a near-infrared wavelength (for example, 800 to 1000 nm).
  • the light source 10 is, for example, a wavelength sweep light source such as a very high frequency pulse laser or an SLD (Superluminescent Diode).
  • the light emitted from the light source 10 is guided inside the optical fiber F and incident on the beam splitter 20.
  • the beam splitter 20 separates the incident light into light that is guided toward the collimator 30a (hereinafter, reference light) and light that is guided toward the collimator 50a (hereinafter, measurement light).
  • the beam splitter 20 is, for example, a cube beam splitter.
  • the collimator 30a changes the reference light guided from the beam splitter 20 into parallel light, and emits the parallel light toward the collimator 30b.
  • the collimator 30b collects the parallel light incident from the collimator 30a and emits the collected reference light toward the reference mirror 40.
  • the collimator 30b incidents the reference light reflected by the reference mirror 40, changes it into parallel light, and emits the parallel light toward the collimator 30a.
  • the collimator 30a collects the parallel light incident from the collimator 30b and guides it toward the beam splitter 20.
  • the collimator 50a changes the measurement light guided from the beam splitter 20 into parallel light, and emits the parallel light toward the galvanometer mirror 60a.
  • the parallel light incident from the collimator 50a is reflected and emitted toward the collimator 50b, respectively.
  • the collimator 50b collects the parallel light incident from the collimator 50a via the galvanometer mirrors 60a and 60b, and irradiates the sample Sm with the collected measurement light.
  • the measurement light applied to the sample Sm is reflected by the reflecting surface of the sample Sm and incident on the collimator 50b.
  • the reflecting surface is not limited to, for example, the boundary surface between the sample Sm and the surrounding environment (for example, the atmosphere) of the sample Sm, and can be a boundary surface between materials or structures having different refractive indexes inside the sample Sm.
  • the light reflected on the reflecting surface of the sample Sm and incident on the collimator 50b is referred to as reflected light.
  • the collimator 50b emits the incident reflected light toward the galvanometer mirror 60b. It is reflected on the surfaces of the galvanometer mirrors 60b and 60a, respectively, and emitted toward the collimator 50a.
  • the collimator 50a collects the parallel light incident from the collimator 50a via the galvanometer mirrors 60a and 60b, and guides the collected reflected light toward the beam splitter 20.
  • the beam splitter 20 guides the reference light reflected by the reference mirror 40 and the reflected light reflected by the sample Sm to the spectroscope 70 via the optical fiber F.
  • the spectroscope 70 includes a diffraction grating and a light receiving element inside the spectroscope 70.
  • the diffraction grating disperses the reference light and the reflected light guided from the beam splitter 20.
  • the separated reference light and reflected light interfere with each other and become interference light that interferes with each other.
  • the light receiving element detects the interference light and generates a signal based on the detected interference light (hereinafter referred to as a detection signal).
  • the light receiving element is configured by arranging a plurality of pixels two-dimensionally on a predetermined imaging surface. The light receiving element outputs the generated detection signal to the signal processing device 100.
  • the signal processing device 100 acquires an OCT signal representing the state of the sample Sm based on the detection signal input from the spectroscope 70.
  • the signal processing device 100 sequentially changes the observation points in a predetermined order, and sequentially obtains detection signals at the individual observation points or for each of a plurality of observation points forming a predetermined unit among the observation points. Cumulatively generate OCT signals within a predetermined area.
  • the signal processing device 100 converts the corrected signal value for each sample point into a luminance value for each pixel, and generates image data indicating the luminance value for each pixel corresponding to each sample point. In the present embodiment, the signal processing device 100 performs the following processing on the acquired OCT signal before generating image data.
  • the signal processing device 100 calculates the phase gradient of the OCT signal for each sample point arranged on a predetermined plane.
  • the signal processing device 100 defines a sample point whose bulk phase error has been determined as a starting point, a sample point whose bulk phase error has not been determined as an end point, and sets a plurality of paths from the starting point to one end point. Then, the signal processing device 100 integrates the phase gradient of each sample point on the set path along each path from the start point to the end point to calculate a bulk phase error (hereinafter, path-specific phase error). To do.
  • the signal processing device 100 synthesizes the path-specific phase errors between the set paths and determines the bulk phase error at the end point.
  • the signal processing device 100 executes a process of determining a bulk phase error for each of the sample points as end points in the plane.
  • the bulk phase error means an error included in the phase of the OCT signal due to the movement of the sample or the fluctuation of the atmosphere or temperature in which the light propagates.
  • the signal processing device 100 corrects the phase of the signal value by subtracting the bulk phase error from the phase of the signal value for each sample point indicating the OCT signal.
  • FIG. 2 is a block diagram showing a functional configuration example of the signal processing device 100 according to the present embodiment.
  • the signal processing device 100 includes a control unit 110 and a storage unit 180.
  • the control unit 110 reads, for example, a program stored in the storage unit 180 in advance by a processor such as a CPU (Central Processing Unit), performs processing instructed by a command described in the read program, and performs its function. Play.
  • a processor such as a CPU (Central Processing Unit)
  • performing a process instructed by a command described in a program may be referred to as executing a program, executing a program, or the like.
  • a part or all of the control unit 110 may be composed of hardware such as LSI (Large Scale Integration), ASIC (Application Specific Integrated Circuit), or dedicated hardware.
  • the control unit 110 includes an optical system control unit 120, a detection signal acquisition unit 130, a phase gradient calculation unit 140, a bulk phase error calculation unit 150, a signal correction unit 160, and an image processing unit 170.
  • the optical system control unit 120 drives a drive mechanism that changes the positions of the galvanometer mirrors 60a and 60b, and scans the sample point which is the observation point of the sample Sm.
  • the sample points of the sample Sm are scanned in each of the depth direction of the sample Sm and the direction of intersection (for example, the vertical direction) in the direction intersecting the depth direction.
  • the depth direction of the sample may be referred to as the z direction in the three-dimensional Cartesian coordinate system, and the first and second directions orthogonal to the z direction may be referred to as the x direction and the y direction, respectively.
  • Signal acquisition in the depth direction of the sample is called A-scan
  • signal acquisition in the direction intersecting that direction (for example, x-direction or y-direction) is called B-scan (B-scan).
  • scan scan
  • the detection signal acquisition unit 130 sequentially acquires detection signals from the spectroscope 70.
  • the detection signal acquisition unit 130 uses the distribution of the reflected light (reflected light) of the measurement light in the depth direction of the sample Sm as a signal value based on the acquired detection signal as a signal value at a predetermined interval in the three-dimensional space where the sample Sm exists. Calculate for each sample point arranged in.
  • the depth direction of the sample Sm corresponds to the incident direction (z direction) of the measurement light.
  • the calculated refractive index distribution corresponds to the intensity distribution of the reflected light.
  • the detection signal acquisition unit 130 repeats the process of calculating the intensity distribution of the reflected light for each observation point changed by signal acquisition along the plane intersecting in the z direction (for example, the xy plane).
  • the detection signal acquisition unit 130 can acquire data (hereinafter, OCT signal) representing the distribution of reflected light in the observable three-dimensional region (hereinafter, observable region).
  • OCT signal indicates a signal value for each sample point distributed at regular intervals in the observable region.
  • the OCT signal is three-dimensional data representing the state of the sample Sm in the observable region, and is also called an OCT volume, a three-dimensional OCT signal (FIG. 3), or the like.
  • the detection signal acquisition unit 130 stores the acquired OCT signal in the storage unit 180.
  • the planes intersecting in the z direction do not necessarily have to be orthogonal to the z direction and may not be parallel to the z direction.
  • the individual sample points do not necessarily have to be arranged on each lattice point of the orthogonal lattice.
  • the individual sample points may be arranged, for example, on each grid point of the diagonal grid or may be arranged aperiodically.
  • the distance between adjacent sample points shall be less than or equal to the spatial resolution of the optical system. In other words, there is no significant difference in the phase component (sample phase) due to the structure of the sample between adjacent sample points.
  • the phase component is canceled and the bulk phase error component is left.
  • signal acquisition involving a change in the target sample point or batch signal acquisition for a plurality of sample points is referred to as "scan", but it is not necessarily a member constituting the optical system (for example, a collimator). It may also include cases without scanning driven by 50b, sample Sm) or detector.
  • the detection signal acquisition unit 130 is responsible for the function related to signal acquisition, and is omitted in the optical system control unit 120. For example, when performing an A-scan at a certain point on the xy plane using FD-OCT as an OCT method, the detection signal acquisition unit 130 converts reflected light having different frequency components depending on the depth of the sample Sm. Based on the interference light, it is detected as a detection signal.
  • the detection signal acquisition unit 130 can Fourier transform the detection signal to calculate the conversion coefficient for each frequency, and determine the conversion coefficient of the depth associated with the frequency as the signal value at that depth. Further, when scanning by driving in the x-direction or the y-direction is not required, the optical system control unit 120 may be omitted.
  • the phase gradient calculation unit 140 reads the OCT signal from the storage unit 180, performs phase unwrapping on the signal value for each sample point indicated by the read OCT signal, and performs phase unwrapping on the signal value.
  • the phase is adjusted to calculate the adjusted phase (unwrapped phase).
  • Phase unwrapping is a process of releasing the limitation so as to ensure spatial continuity for the phase limited to the range of ⁇ to ⁇ .
  • the phase gradient calculation unit 140 calculates the phase gradient for each sample point arranged on a predetermined two-dimensional plane based on the adjusted phase for each sample point (FIG. 3: phase gradient analysis S01).
  • the phase gradient is a two-dimensional vector including partial derivatives with respect to each of the first direction (for example, the x direction) and the second direction (for example, the y direction) that are parallel to the two-dimensional plane and orthogonal to each other.
  • the two-dimensional plane may be a plane orthogonal to the z direction or a plane not orthogonal to each other as long as the planes intersect in the z direction. Further, the two-dimensional plane may be a plane that crosses the observable area as long as it can be projected from an arbitrary point in the three-dimensional observable area, or may be a virtual plane. ..
  • phase gradient calculation unit 140 stores in the storage unit 180 data indicating the phase gradient calculated for each sample point in the two-dimensional region as phase gradient data (FIG. 3: two-dimensional phase gradient Pg01).
  • the bulk phase error calculation unit 150 reads the phase gradient data from the storage unit 180, performs the processing described below based on the phase gradient for each sample point indicated by the read phase gradient data, and arranges them in the two-dimensional region. Calculate the bulk phase error at each sample point.
  • the bulk phase error calculation unit 150 starts from one or a plurality of sample points for which the bulk phase error has been determined, and extends from the start point to the end point where the bulk phase error is one of the undetermined sample points. Set the route.
  • the bulk phase error calculation unit 150 integrates the phase gradient for each sample point for each of the set paths from the start point to the end point along each path (Fig. 3: Phase estimation S02 by path integral) and integrates at the end point. Calculate the value as a phase error for each path.
  • the bulk phase error calculation unit 150 synthesizes the combined value calculated by synthesizing the phase errors for each path calculated for each path between a plurality of paths set, and the bulk phase error at the end point (FIG. 3: bulk phase error Bp01). Determined as. However, the bulk phase error calculation unit 150 sets one of the sample points on the two-dimensional plane as a reference point (for example, the origin of the two-dimensional region), and sets the bulk phase error at the reference point as a predetermined reference value (for example, for example). , 0) in advance. The bulk phase error calculation unit 150 stores the bulk phase error data indicating the bulk phase error for each defined sample point in the storage unit 180.
  • the signal correction unit 160 reads the bulk phase error data and the OCT signal from the storage unit 180.
  • the signal correction unit 160 performs a process of correcting the phase of the signal value for each sample point indicated by the OCT signal based on the bulk phase error for each sample point indicated by the read bulk phase error data.
  • the signal correction unit 160 is a sample in the observable region in which the positions of the sample points (hereinafter, two-dimensional sample points) on the two-dimensional region indicated by the bulk phase error data and the positions of the points projected on the two-dimensional plane are common.
  • a point hereinafter referred to as a three-dimensional sample point
  • the number of three-dimensional sample points specified can be plural for one two-dimensional sample point.
  • the signal correction unit 160 corrects the phase of the signal value at the specified three-dimensional sample point by subtracting the bulk phase error of the corresponding two-dimensional sample point.
  • the signal correction unit 160 may further perform known correction processing (for example, refocusing, speckle suppression, etc.) on the signal value whose phase is corrected for each three-dimensional sample point. Good.
  • the signal correction unit 160 stores in the storage unit 180 a correction OCT signal indicating a signal value for each three-dimensional sample point obtained by the correction process.
  • the image processing unit 170 reads the corrected OCT signal from the storage unit 180 and converts the signal value for each three-dimensional sample point into a luminance value using a predetermined conversion function.
  • the converted luminance value takes a value within a value range that can be expressed by the bit depth of each pixel.
  • Image data having a luminance value for each sample point on a two-dimensional plane to be observed (hereinafter referred to as an observation plane) as a luminance value for each pixel is generated.
  • sample points may be referred to as pixels.
  • the image processing unit 170 may define a cross section that crosses at least a part of the observable region as a two-dimensional plane to be observed based on an operation signal input from an operation input unit (not shown).
  • the operation signal indicates, for example, the depth from the sample surface, the observation direction with respect to the sample surface, and the like as the requirements of the observation target plane.
  • the operation input unit may include, for example, a button, a knob, a dial, a mouse, a joystick, and other members that accept user operations and generate operation signals according to the received operations.
  • the operation input unit may be an input interface that receives an operation signal wirelessly or by wire from another device (for example, a portable device such as a remote controller).
  • the image processing unit 170 may output the generated image data to a display unit (not shown) and display an OCT image based on the image data.
  • the display unit may be, for example, a liquid crystal display (LCD), an organic electroluminescence display (OLED), or the like.
  • the image processing unit 170 may store the generated image data in the storage unit 180, or may wirelessly or wire-transmit it to another device via a communication unit (not shown).
  • the storage unit 180 stores various data used for processing executed by the control unit 110 and various data acquired by the control unit 110.
  • the storage unit 180 includes, for example, a non-volatile (non-temporary) storage medium such as a ROM (Read Only Memory), a flash memory, and an HDD (Hard Disk Drive).
  • the storage unit 180 includes, for example, a volatile storage medium such as a RAM (Random Access Memory) and a register.
  • the phase gradient calculation unit 140 is a two-dimensional sample point (x, y, z) arranged in a two-dimensional region based on the adjusted phase ⁇ (x, y, z) for each three-dimensional sample point (x, y, z).
  • the phase gradient (Bx, By) is calculated for each y).
  • the three-dimensional sample points (x, y, z) are distributed in the observable region of the rectangular parallelepiped at predetermined intervals in the x direction, the y direction, and the z direction, respectively.
  • the x-direction, y-direction, and z-direction are orthogonal to each other and correspond to the directions of the respective sides forming the outer edge of the observable region.
  • the two-dimensional region is a two-dimensional xy plane parallel to the x and y directions.
  • the coordinates of the two-dimensional sample points (x, y) that project the three-dimensional sample points (x, y, z) onto the two-dimensional region are the coordinates of the three-dimensional sample points (x, y, z) in the z direction. It can be obtained by ignoring the value.
  • the coordinates of the three-dimensional sample point (x, y, z) corresponding to the two-dimensional sample point (x, y) are the x-coordinate value and the y-coordinate value of the two-dimensional sample point (x, y). These are sample points indicated by the common x-coordinate value and y-coordinate value, respectively.
  • the phase gradient calculation unit 140 is, for example, a three-dimensional sample adjacent to the spectrum A (x, y, z) of the signal value at the three-dimensional sample point (x, y, z) in the x direction.
  • the argument of the sum of the values multiplied by the complex conjugate A * (x + ⁇ x, y, z) of the spectrum of the signal value at the point (x + ⁇ x, y, z) over the z direction is the corresponding two-dimensional sample point ( It is calculated as the x-direction component B x (x, y) of the phase gradient at x, y).
  • Phase gradient calculating unit 140 the corresponding two-dimensional sample points (x, y) y direction component B y (x, y) of the phase gradient in regard to the three-dimensional sample points as shown in equation (2) (x, Complex conjugate A * (x, y + ⁇ y, z) of the signal value spectrum at the signal value spectrum A (x, y, z) at y, z) and the signal value spectrum at the three-dimensional sample point (x, y + ⁇ y, z) adjacent in the y direction. ) Can be calculated.
  • ⁇ x and ⁇ y indicate the intervals of the three-dimensional sample points in the x and y directions, respectively, and indicate the intervals of the two-dimensional sample points.
  • arg indicates the argument of the complex number ...
  • phase gradient B calculated for each two-dimensional sample point is used for calculating the bulk phase error as described above.
  • the integral path is not necessarily limited to the portion parallel to the x direction or the y direction, but may include a portion having a direction intersecting both the x direction and the y direction.
  • the phase gradient calculation unit 140 sets the next sample point (x + ⁇ x, y + ⁇ y) and the phase gradient B xy (x, y) at the sample points (x, y) of the portion based on the equation (3). ) Can be calculated.
  • the bulk phase error calculation unit 150 presets a plurality of integral paths having individual phase determination target pixels as end points for each two-dimensional sample point (hereinafter, also referred to as phase determination target pixel) to be determined for the phase error. I will do it.
  • the starting point of each integral path is a two-dimensional sample point whose phase error has been determined (hereinafter, also referred to as a phase-determined pixel). However, the starting point is a two-dimensional sample point located within a predetermined range from the ending point. This limits the number of integral paths and does not become excessive.
  • the bulk phase error calculation unit 150 integrates the phase gradient B for each sample point (x, y) on the path Pt from the start point Op to the end point Dp of each path Pt, for example, as shown in the equation (4).
  • the integral value at the end point obtained by this is calculated as the phase error ⁇ e and Pt for each path.
  • B (x, y) represents the phase gradient, the sample points (x, y) x-direction component B x (x, y) in the y direction component B y (x, y) elements It is a two-dimensional vector including as. n (x, y) indicates a direction vector in the tangential direction of the path Pt. More specifically, as a direction vector n (x, y), a vector having a sample point (x, y) as a starting point and an adjacent sample point on the path Pt as an ending point can be used. • indicates the inner product of the vectors. Then, the bulk phase error calculation unit 150 synthesizes the path-specific phase errors ⁇ e and Pt calculated for the path Pt between the paths having a common end point Dp, and calculates the bulk phase error ⁇ e at the end point Dp.
  • FIG. 4 is an explanatory diagram showing an example of the integral path according to the present embodiment.
  • four paths Pt01-Pt04 are set for the phase-determined pixel Op00 as the starting point and the phase-determining target pixel Dp11 as the ending point.
  • the coordinates of the phase-determined pixel Op00 and the phase-determined pixel Dp11 on the two-dimensional plane are (0,0) and (1,1), respectively.
  • the end point is adjacent to the start point in an oblique direction that intersects both the x-coordinate and the y-coordinate.
  • the x-coordinate and y-coordinate are normalized so that the intervals ⁇ x and ⁇ y are 1, respectively.
  • the paths Pt01-Pt04 are formed by sequentially connecting elementary pieces starting from one sample point and ending at adjacent sample points in the x-direction or the y-direction.
  • the vector showing each element is used as the direction vector shown in the equation (4) for calculating the phase error for each path.
  • the bulk phase error calculation unit 150 determines the number of pieces forming each path Pt01-Pt04 as the path length (hereinafter, the path length).
  • the path length the path length
  • the path lengths of the paths Pt01, Pt02, Pt03, and Pt04 are 2, 2, 4, and 4, respectively.
  • the bulk phase error calculation unit 150 determines the path-specific phase errors ⁇ e, Pt01 , ⁇ e, Pt02 , ⁇ e, Pt03 , ⁇ e, Pt 04 related to the paths Pt01, Pt02, Pt03, and Pt04 in the phase determination target pixel Dp11.
  • each B x (0,0) + B y (1,0), B y (0,0) + B x (0,1), B x (0,0) + B x (1,0) + B y (2 , 0) -B x (2,1) , B y (0,0) + B y (0,1) + B x (0,2) can be calculated as -B y (1,2).
  • the bulk phase error calculation unit 150 calculates, for example, a simple average value between paths as a bulk phase error in which path-specific phase errors are combined between paths.
  • a process of dividing the total sum of the phase errors for each path by 4 which is the number of paths and normalizing the phase error is included.
  • the bulk phase error calculation unit 150 may calculate the weighted average value of the phase error for each path as the bulk phase error by using a value depending on the path length as the weighting coefficient for each path. For example, the bulk phase error calculation unit 150 normalizes the weighting coefficient for each path to ap (a is a real number larger than 0 and less than 1 and p is the path length), and the sum of the paths is normalized to 1. The real number given can be used. In that case, in the example shown in FIG.
  • the path-specific phase error ⁇ e, Pt01, ⁇ e, Pt02, ⁇ e, Pt03, weight factor for ⁇ e, Pt04 respectively a 2/2 (a 2 + a 4), a 2/2 (a 2 + a 4), a 4/2 (a 2 + a 4), the a 4/2 (a 2 + a 4).
  • the bulk phase error calculation unit 150 may use a value that is inversely proportional to the path length p and the sum of the paths is normalized to 1 as the weighting coefficient for each path. In that case, in the example shown in FIG.
  • the weighting coefficients for the path-specific phase errors ⁇ e, Pt01 , ⁇ e, Pt02 , ⁇ e, Pt03 , ⁇ e, and Pt 04 are 1/3, 1/3, and 1/6, respectively. , 1/6. According to these methods, it is possible to reduce the contribution of the path-specific phase error of the path having a large path length to the bulk phase error.
  • the bulk phase error calculation unit 150 may calculate the bulk phase error ⁇ b (t) (x 0 , y 0 ) at the end point (x 0 , y 0) using the equation (5).
  • ⁇ e, Pt (x 0 , y 0 ) indicates a path-specific phase error with respect to the path Pt.
  • the sum vector which is the sum of the unit vectors having a magnitude of 1 and having the declination of the phase error for each path, is calculated, and the declination of the calculated sum vector is the bulk phase error ⁇ b.
  • T Indicates that it is defined as (x 0 , y 0).
  • various integration paths can be set according to the positional relationship between the phase-determined pixel and the phase-determined target pixel.
  • a set of a plurality of integral paths corresponding to the positional relationship between the phase-determined pixel and the phase-determined target pixel is referred to as a path pattern.
  • the route pattern illustrated in FIG. 4 is called a route pattern Pn00.
  • 5A-5I are diagrams showing other examples of the route pattern according to the present embodiment.
  • the end point is the phase determination target pixel Dp11
  • the position of the phase determined pixel as the start point is x with respect to the phase determination target pixel Dp11.
  • the number of phase-determined pixels as a starting point is not limited to one, and may be two or more.
  • the number of phase-determined pixels as the starting point is one.
  • the number of phase-determined pixels as the starting point is two.
  • the number of phase-determined pixels as the starting point is three.
  • the path pattern Pn01 illustrated in FIG. 5A has three paths Pt11-Pt13. Each of the paths Pt11 to Pt13 starts from the phase-determined pixel Op01. The path lengths of the paths Pt11, Pt12, and Pt13 are 3, 1, and 3, respectively.
  • the path pattern Pn02 illustrated in FIG. 5B has three paths Pt21-Pt23. The paths Pt21 to Pt23 each start from the phase-determined pixel Op10. The path lengths of the paths Pt21, Pt22, and Pt23 are 3, 1, and 3, respectively.
  • the path pattern Pn03 illustrated in FIG. 5C has two paths Pt31 and Pt32. The paths Pt31 and Pt32 start from the phase-determined pixels Op01 and Op02, respectively. The path lengths of the paths Pt31 and Pt32 are 1, 3 respectively.
  • the path pattern Pn04 illustrated in FIG. 5D has two paths Pt41 and Pt42.
  • the paths Pt41 and Pt42 start from the phase-determined pixels Op10 and Op20, respectively.
  • the path lengths of the paths Pt41 and Pt42 are 1, 3 respectively.
  • the path pattern Pn05 illustrated in FIG. 5E has four paths Pt51-Pt54.
  • the paths Pt51, Pt52, Pt53, and Pt54 start from the phase-determined pixels Op01, Op10, Op01, and Op10, respectively.
  • the path lengths of the paths Pt51, Pt52, Pt53, and Pt54 are 1, 1, 3, and 3, respectively.
  • the path pattern Pn06 illustrated in FIG. 5F has three paths Pt61-Pt63.
  • the paths Pt61, Pt62, and Pt63 start from the phase-determined pixels Op00, Op01, and Op02, respectively.
  • the path lengths of the paths Pt61, Pt62, and Pt63 are 2, 1,
  • the path pattern Pn07 illustrated in FIG. 5G has three paths Pt71-Pt73.
  • the paths Pt71, Pt72, and Pt73 start from the phase-determined pixels Op00, Op10, and Op20, respectively.
  • the path lengths of the paths Pt71, Pt72, and Pt73 are 2, 1, and 2, respectively.
  • the path pattern Pn08 illustrated in FIG. 5H has four paths Pt81-Pt84.
  • the paths Pt81, Pt82, Pt83, and Pt84 start from the phase-determined pixels Op01, Op10, Op01, and Op20, respectively.
  • the path lengths of the paths Pt81, Pt82, Pt83, and Pt84 are 1, 1, 3, and 2, respectively.
  • 5I has four paths Pt91-Pt94.
  • the paths Pt91, Pt92, Pt93, and Pt94 start from the phase-determined pixels Op01, Op02, Op10, and Op10, respectively.
  • the path lengths of the paths Pt91, Pt92, Pt93, and Pt94 are 1, 3, 1, and 3, respectively.
  • a route pattern having a route in which the x-coordinate of the end point is less than the x-coordinate of the start point or the y-coordinate of the end point is less than the y-coordinate of the start point may be set, or a route pattern having a larger number of start points may be set. ..
  • 6A-6D and 7A-7C show an example of a set of start and end points of a path pattern according to this embodiment.
  • 6A-6D and 7A-7C show all sets that can be set in the region in the xy plane of 3 samples x 3 samples centered on the end point.
  • a set and a set whose origin distribution is symmetric with respect to the x-axis or y-axis passing through the end point and a set whose rotation symmetry is about the end point are regarded as a common set and are not shown.
  • the starting point and the ending point are indicated by diagonal lines and shading, respectively.
  • the illustration of the route is omitted.
  • 6A-6D show a set in which the number of starting points is 1-4, respectively.
  • the number of starting points is 1, 2, 3, or 4, 2 sets, 4 sets, 8 sets, and 12 sets can be set, respectively.
  • 7A-7C show a set in which the number of starting points is 5-7, respectively. When the number of starting points is 5, 6 and 7, 8 sets, 4 sets and 2 sets can be set respectively.
  • the starting point is a phase-determined image adjacent to one or both (diagonal directions) of the x-direction and the y-direction from the ending point, but the starting point is predetermined from the ending point. If it is within the range, it is not limited to this.
  • the starting point may be a phase-determined pixel within 2 samples from the ending point.
  • all of the plurality of integral paths are within a predetermined path setting range (for example, FIG. 8B, FIG. 8C, FIG. 9: path setting range Pa) from the end point. By including it and setting it so that part or all of it does not fall outside the predetermined range, it is possible to prevent the integral path from becoming excessive.
  • the integral path is formed by sequentially connecting elementary pieces (links) consisting of sample points adjacent in the x-direction or y-direction from a certain sample point from the start point to the end point.
  • the starting point of each piece may be a sample point arranged within the spatial resolution of the optical system from the end point of the piece.
  • the element piece may be, for example, a section connecting two diagonally adjacent sample points that intersect both the x-direction and the y-direction, or straddle one or more sample points without being adjacent to each other (). Skip) It may be a section connecting two sample points.
  • sample points in the middle of the integral path may be sample points for which the bulk phase error has already been determined, or sample points for which the bulk phase error has not been determined. It may be.
  • the phase gradient calculation unit 140 determines the phase difference of the OCT signal between each sample point in the two-dimensional region and the end point using the sample point as the start point of the element piece for each position of the end point that can be taken as a piece. Calculate as a phase gradient.
  • a plurality of path patterns may be set in advance in the bulk phase error calculation unit 150.
  • the bulk phase error calculation unit 150 preliminarily sets a path pattern for each sample point in the two-dimensional region, starting from the sample point and satisfying the positional relationship and the correspondence relationship with the sample points for which the bulk phase error has already been determined as the starting point. Select from the set route patterns.
  • the bulk phase error calculation unit 150 calculates the bulk phase error at the end point by the above method with each of the plurality of paths indicated by the selected path patterns as the integral path. Then, the bulk phase error calculation unit 150 changes the target sample point for which the bulk phase error is calculated to another sample point whose bulk phase error has not been determined.
  • the bulk phase error calculation unit 150 may repeat the calculation of the bulk phase error and the change of the target sample points until the bulk phase error is determined for all the sample points in the two-dimensional region.
  • the target sample point corresponds to the above-mentioned phase determination target pixel.
  • the bulk phase error calculation unit 150 is provided with a calculation order (iteration) for calculating the bulk phase error so that one sample point in the two-dimensional region is set as a reference point and the sample points closer to the reference point are prioritized. It may be determined in advance. Further, when the bulk phase error calculation unit 150 uses a sample point whose end point is a sample point whose calculation order is the mth time (m is an integer of 1 or more), the bulk phase error calculation unit 150 is a sample point whose calculation order is the m-1th time. May be used as a starting point, and a route pattern corresponding to the positional relationship between the ending point and the starting point may be selected.
  • the reference point is set as the sample point whose operation order is the 0th time
  • the bulk phase error is set as a predetermined reference value.
  • the number of sample points having the same calculation order is not necessarily one, but may be two or more.
  • the number of sample points having the same calculation order may be two or more, but the calculation order of the bulk phase error may be a predetermined order or randomly determined between two or more sample points in the m-th time. You may.
  • the bulk phase error calculation unit 150 has the origin (0,0) of the two-dimensional region as the reference point Op00, and is sandwiched in the x-direction and the y-direction from the origin. Set the calculation order in the order of the distances. Sample points in the same calculation order are distributed on lines that are orthogonal to each other in the diagonal direction. As shown in FIG. 8B, as the sample points having the first calculation order, the sample points located at (1,0) and (0,1) are set as the target sample points Dp10 and Dp01, respectively. At this time, the bulk phase error calculation unit 150 can select the path patterns Pn01 and Pn02 for the target sample points Dp10 and Dp01, respectively.
  • the bulk phase error calculation unit 150 rejects the selected path patterns Pn01 and Pn02 as invalid paths because the paths Pt11 and Pt21 each pass outside the range of the two-dimensional region in which the sample points are arranged. , The remaining valid paths are used to calculate the bulk phase error.
  • the weighting coefficient used when synthesizing the phase error for each route between effective routes may be predetermined based on the number of effective routes or the route length of the effective route.
  • the sample points whose operation order is the second are the target sample points Dp11, Dp20, and Dp02 located at (1,1), (2,0), and (1,0), respectively. Is set as.
  • the bulk phase error calculation unit 150 can select the path patterns Pn05, Pn01, and Pn02 for each of the target sample points Dp10 and Dp01.
  • the bulk phase error calculation unit 150 rejects the selected path patterns Pn01 and Pn02 as invalid paths because the paths Pt11 and Pt21 each pass outside the range of the two-dimensional region in which the sample points are arranged. , The remaining valid paths are used to calculate the bulk phase error. As shown in FIG.
  • the sample points located in (2, 1) and (1, 2) are set as the target sample points Dp21 and Dp12, respectively.
  • the sample point whose calculation order is the third time includes the target sample point Dp22 farthest from the reference point Op00 in the two-dimensional region.
  • the bulk phase error calculation unit 150 can select the path patterns Pn05, Pn05, and Pn00 for each of the target sample points Dp21, Dp12, and Dp22.
  • sample points are arranged in the paths Pt54, Pt53, Pt03, and Pt04, respectively. Since it passes outside the range of the two-dimensional region, it is rejected as an invalid path, and the remaining valid path is used for calculating the bulk phase error.
  • the calculation order between the sample points in the two-dimensional region is not limited to the examples shown in FIGS. 8A-8C and 9A.
  • the sample point as the reference point may be set as the center of the two-dimensional region, and the calculation order may be set so that the sample around the reference point that is closer to the reference point has a higher calculation order.
  • the bulk phase error calculation unit 150 sets the target sample point whose end point is the sample point whose calculation order is mth (m is an integer of 1 or more)
  • the sample point whose calculation order is m-1th May be used as a starting point, and a route pattern corresponding to the positional relationship between the ending point and the starting point may be selected.
  • FIG. 9B shows, as an example, when calculating the bulk phase error at the sample point A in which the path Sq indicating the calculation order is set in a clockwise spiral shape centered on the reference point.
  • Bulk phase error calculating unit 150 sets a region of 3 ⁇ 3 points around the sample point A as a route setting range P A. Since the phase determination already pixels in the path setting range P A is only a reference point, bulk phase error calculating unit 150 selects the route pattern illustrated in Figure 5B.
  • the bulk phase error calculation unit 150 sets a region of 3 ⁇ 3 points centered on the sample point B as the route setting range P B.
  • the bulk phase error calculation unit 150 is illustrated at the left end of the first row in FIG. 6D. Select the route pattern to be used.
  • the calculation order for each sample is determined so that the closer the distance from the reference point is, the priority is given, and the path integral starting from the immediately preceding calculation order is executed.
  • the number of operations until the phase error is calculated is reduced. Therefore, the accumulation of errors is suppressed in the entire two-dimensional region.
  • the bulk phase error calculation unit 150 may predetermine any one of a plurality of path patterns, or randomly select one of the plurality of path patterns. May be good.
  • the bulk phase error calculation unit 150 plots each phase-determined pixel as a bulk phase error calculation process (hereinafter, line process) in each row until the phase-determined pixel reaches the sample point at the end of the two-dimensional region.
  • the bulk phase error is calculated using the path pattern shown in 5A, and the process of changing the phase-determined pixel to an adjacent sample point in the x-axis direction is repeated.
  • the bulk phase error calculation unit 150 calculates the bulk phase error at the sample point at the start of the next row adjacent in the y-axis direction using the path pattern shown in FIG. 5B. Then, the bulk phase error calculation unit 150 may repeat the line processing until the line to be determined for the bulk phase reaches the end of the two-dimensional region.
  • FIG. 10A-10F are diagrams showing an example of phase correction of the OCT signal according to the present embodiment.
  • FIG. 10A shows an example of an OCT image based on an OCT signal.
  • the OCT image shown in FIG. 10A shows the state of the observation target plane forming a part of the observable region.
  • This OCT image is based on a raw OCT signal for which bulk phase errors have not yet been compensated.
  • FIG. 10D shows a spectrum obtained by converting the OCT signal according to FIG. 10A into a spatial frequency domain.
  • the vertical axis and the horizontal axis indicate the spatial frequencies in the x and y directions, respectively, and the center indicates the origin.
  • the intensity of the spectrum is represented by shading. The brighter the part, the higher the intensity, and the darker the part, the lower the intensity.
  • the intensity is asymmetric with respect to the origin, and the center of the frequency band having the highest intensity is displaced from the origin in the x direction. In addition, the spread of the high-strength portion is relatively large.
  • FIG. 10B shows an example of a bulk phase error calculated using the conventional phase difference method.
  • the phase is represented by shading. The brighter the part, the closer the phase is to ⁇ , and the darker the part, the closer the phase is to - ⁇ .
  • the bulk phase error changes relatively regularly with respect to the x direction, but changes in the y direction are more irregular than in the x direction.
  • FIG. 10E shows a spectrum corrected using the bulk phase error shown in FIG. 10B.
  • the intensity of the spectrum is symmetric with respect to the y-axis passing through the origin in the x-direction, and the spread in the x-direction is relaxed. However, the spread of the spectral intensity in the y direction is rather widened. This indicates that the correction in the y direction is inappropriate.
  • FIG. 10C shows an example of a bulk phase error calculated by the signal processing apparatus 100 according to the present embodiment.
  • the phase is represented by shading.
  • This bulk phase error was calculated using the method described with reference to FIGS. 8A-8C and 9.
  • the bulk phase error changes regularly in both the x-direction and the y-direction, and shows a relatively smooth striped pattern.
  • FIG. 10F shows the spectrum of the OCT signal corrected using the bulk phase error shown in FIG. 10C.
  • the signal correction unit 160 can remove the bulk phase error included in the OCT signal by correcting the OCT signal using the equation (6).
  • Equation (6) is obtained by subtracting the bulk phase error ⁇ b (t) (x, y) from the phase of the spectrum A (x, y, z) of the signal value at the three-dimensional sample point (x, y, z). It is shown that the spectrum A'(x, y, z) of the corrected signal value is calculated. The intensity of the spectrum is symmetric with respect to the origin in both the x-direction and the y-direction, and the spread is relaxed. This indicates that the present embodiment can alleviate the directional dependence caused by the phase error as compared with the conventional phase difference method, and can realize more appropriate spectral correction.
  • FIG. 11A and 11B show that the average value of the original data varies significantly among the samples, and that the average value of the simple integral and the proposed method is almost 0, which is significantly smaller than that of the original data.
  • the average value of the spatial frequency f y, towards simple integration than the proposed method is rather variations are suppressed.
  • FIG. 11D standard deviation of the spatial frequency f y is increased from the original data by simple integration. This indicates that the spectrum spreads in the y direction rather than the x direction.
  • the proposed method no significant change in the standard deviations of the spatial frequencies f x and f y is observed with the original data. This indicates that according to the proposed method, the bulk phase error can be corrected without causing a spectral bias in a specific direction.
  • the signal correction unit 160 may refocus the corrected OCT signal (digital refocusing, digital refocusing). In the refocusing, the signal correction unit 160 executes the processing of each of the following steps.
  • the signal correction unit 160 performs a two-dimensional discrete Fourier transform (DFT) on the corrected OCT signal (for example, front OCT signal, en face OCT slice) in the two-dimensional region. Calculate the complex spatial frequency spectrum in the spatial frequency domain.
  • Step S12 the signal correction unit 160, the spatial frequency (f x, f y) constituting the complex spectral transform coefficients for each phase filter H -1 (f x, f y) for the spatial frequency multiplying the Calculates the filtered spectrum including the multiplication value obtained by multiplication.
  • Phase filter H -1 (f x, f y ) is an example of, shown in Equation (7).
  • f r represents the radial direction of the spatial frequency.
  • f r corresponds to the square root of the sum of squares of the space f x in the x direction and the space frequency f y in the y direction.
  • ⁇ c indicates the central wavelength of the incident light.
  • z 0 indicates the amount of defocus. The amount of defocus z 0 generally depends on the optical system.
  • Phase filter H -1 (f x, f y ) is to compensate for the phase shifts conspicuous higher spatial frequency, making it possible to alleviate or eliminate defocus. Such phase shifts typically appear in the OCT image as "blurring".
  • Step S13 The signal correction unit 160 performs a two-dimensional inverse discrete Fourier transform (IDFT) on the filtered spectrum to generate the filtered OCT signal as a refocused OCT signal. .. Then, the signal correction unit 160 outputs the refocused OCT signal to the image processing unit 170. According to the image processing unit 170, a clearer OCT image can be acquired based on the OCT signal input from the signal correction unit 160.
  • IDFT two-dimensional inverse discrete Fourier transform
  • the signal processing device 100 includes a phase gradient calculation unit 140 and a bulk phase error calculation unit 150.
  • the phase gradient calculation unit 140 calculates the phase gradient of the OCT signal representing the state of the sample for each sample point arranged in a predetermined plane (for example, a two-dimensional region).
  • the bulk phase error calculation unit 150 sets the phase gradient for each sample point for each of the plurality of paths from the starting point where the bulk phase error has been determined to the end point where the bulk phase error has not been determined.
  • the phase error for each path at the end point is calculated by integrating along the path of, and the phase error for each path is combined between a plurality of paths to determine the bulk phase error at the end point.
  • the bulk phase error calculation unit 150 calculates as the bulk phase error, for example, the average value between the paths of the phase error for each path as the bulk phase error. According to this configuration, random components of errors, which may differ depending on the route, are relatively reduced by being synthesized between the routes. Therefore, the non-path component is estimated as a bulk phase error more accurately than simply accumulating the phase gradients.
  • the bulk phase error calculation unit 150 may calculate the weighted average value of the phase error for each path using the weighting coefficient based on the path length of each of the plurality of paths as the average value between the paths. ..
  • the bulk phase error can be calculated more accurately in consideration of the influence of the error depending on the path length. For example, by reducing the weighting coefficient for a route having a longer route length, the contribution of the phase error for each route of the route can be reduced when the influence of the error is stronger as the route length is larger. On the contrary, by increasing the weighting coefficient as the route length is longer, when the influence of the error is weaker as the route length is larger, the contribution of the phase error for each route of the route can be increased.
  • each of the plurality of paths used for calculating the bulk phase error may be included within a predetermined range from the end point. According to this configuration, since each of the plurality of paths is limited to a predetermined range, it is possible to avoid an excessive number of paths related to the calculation of the bulk phase error.
  • each of the plurality of paths used for calculating the bulk phase error is a section from the nth (n is an integer of 1 or more and N-1 or less) sample point to the n + 1th sample point.
  • the nth element (N is an integer of 2 or more) is included, and the n + 1 sample point may be a sample point within the range of the spatial resolution or less of the optical system that acquired the OCT signal from the nth sample point.
  • each of the plurality of paths is formed by connecting one or a plurality of elementary pieces consisting of two sample points that do not cause a significant difference in the signal value of the OCT signal, and therefore, from the calculated bulk phase error. The contribution of the signal component of the OCT signal can be suppressed.
  • a plurality of path patterns indicating a plurality of paths are set in advance in the bulk phase error calculation unit 150 in association with the positional relationship between the start point and the end point, and the bulk phase error is calculated.
  • a plurality of paths represented by a path pattern that satisfies the positional relationship starting from the sample point for which the bulk phase error has been determined may be selected with the target sample point as the end point.
  • the bulk phase error calculation unit 150 uses one sample point in the two-dimensional plane (for example, a two-dimensional region) as a reference point, and the sample point closer to the reference point is given priority.
  • the order for calculating the bulk phase error is predetermined and the sample points whose specified order is the mth time (m is an integer of 1 or more) are set as the target sample points, the determined order is the m-1th time.
  • a plurality of routes represented by a route pattern satisfying the positional relationship starting from the sample point of the above may be selected.
  • the influence of the accumulated error on the sequentially calculated bulk phase error is suppressed as compared with the case where the operation order is arbitrary. Therefore, the bulk phase error should be calculated with higher accuracy for the entire two-dimensional plane. Can be done.
  • the signal processing device 100 may include a signal correction unit 160.
  • the signal correction unit 160 may subtract the bulk phase error at the end point from the phase of the signal value of the OCT signal at the sample points arranged in the three-dimensional space where the end point and the position projected on the plane are common. As a result, the bulk phase error included in the three-dimensional OCT signal is compensated, so that the spectrum bias and spread due to the bulk phase error can be reduced, and a higher quality optical interference tomographic image can be obtained.
  • the signal processing device 100 is a part of the optical interference tomogram 1 is taken as an example, but the present invention is not limited to this.
  • the signal processing device 100 may be a single device that is independent of the optical interference tomometer 1 and does not have an optical system.
  • the optical system control unit 120 may be omitted in the control unit 110 of the signal processing device 100.
  • the detection signal acquisition unit 130 is not limited to the optical system, and may acquire a detection signal or an OCT signal from another device such as a data storage device or a PC by wire or wirelessly, for example, via a network.
  • the signal processing device 100 may include the above-mentioned operation input unit and display unit, or one or both of them may be omitted.
  • the image processing unit 170 may be omitted in the control unit 110 of the signal processing device 100.
  • the control unit 110 of the signal processing device 100 may output the correction OCT signal generated by the signal correction unit 160 to another device such as a data storage device, a PC, or an image processing device.
  • the output destination device has the same function as the image processing unit 170, that is, a function of generating image data based on the corrected OCT signal input from the signal processing device 100 and displaying an OCT image based on the generated image data. May have.
  • the case where the two-dimensional region in which the sample points related to the calculation of the phase gradient or the bulk phase error are arranged is set on the xy plane is used as an example, but the present invention is not limited to this. It suffices to be set in a two-dimensional plane capable of forming a correspondence between the three-dimensional sample points arranged in the three-dimensional space forming the observable region and the two-dimensional sample points obtained by projecting the three-dimensional sample points.
  • a two-dimensional plane may be, for example, an yz plane, a zx plane, or the like.
  • the three-dimensional sample points (x, y, z) are associated with the two-dimensional sample points (y, z). Further, as the two-dimensional region, an observation target plane in which sample points corresponding to the pixels constituting the OCT image are arranged may be used.
  • the signal correction unit 160 may perform phase correction of the OCT signal based on the bulk phase error for each observation target plane.
  • a part of the signal processing device 100 in the above-described embodiment for example, all or part of the control unit 110 may be realized by a computer.
  • the program for realizing this control function may be recorded on a computer-readable recording medium, and the program recorded on the recording medium may be read by the computer system and executed.
  • the "computer system” referred to here is a computer system built in the signal processing device 100, and includes hardware such as an OS and peripheral devices.
  • the "computer-readable recording medium” refers to a portable medium such as a flexible disk, a magneto-optical disk, a ROM, or a CD-ROM, or a storage device such as a hard disk built in a computer system.
  • a "computer-readable recording medium” is a medium that dynamically holds a program for a short period of time, such as a communication line when a program is transmitted via a network such as the Internet or a communication line such as a telephone line.
  • a program may be held for a certain period of time, such as a volatile memory inside a computer system serving as a server or a client.
  • the above-mentioned program may be a program for realizing a part of the above-mentioned functions, and may be a program for realizing the above-mentioned functions in combination with a program already recorded in the computer system.
  • a part or all of the signal processing device 100 in the above-described embodiment and modification may be realized as an integrated circuit such as an LSI (Large Scale Integration).
  • LSI Large Scale Integration
  • Each functional block of the signal processing device 100 may be individually made into a processor, or a part or all of them may be integrated into a processor.
  • the method of making an integrated circuit is not limited to LSI, and may be realized by a dedicated circuit or a general-purpose processor. Further, when an integrated circuit technology that replaces an LSI appears due to advances in semiconductor technology, an integrated circuit based on this technology may be used.
  • Optical interference tomography 10 ... Light source, 20 ... Beam splitter, 30a, 30b, 50a, 50b ... Collimator, 40 ... Reference mirror, 60a, 60b ... Galvano mirror, 70 ... Spectrometer, 100 ... Signal processor, 120 ... Optical system control unit, 130 ... Detection signal acquisition unit, 140 ... Phase gradient calculation unit, 150 ... Bulk phase error calculation unit, 160 ... Signal correction unit, 170 ... Image processing unit

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

位相勾配算出部は試料の状態を表す光干渉断層信号の光の照射方向に交差する平面における位相勾配を前記平面に配列されたサンプル点ごとに算出し、バルク位相エラー算出部はバルク位相エラーが決定済のサンプル点である起点からバルク位相エラーが未決定のサンプル点である終点までの複数の経路のそれぞれについて、サンプル点ごとの前記位相勾配をそれぞれの経路に沿って積分して前記終点における経路別位相エラーを算出し、前記経路別位相エラーを前記複数の経路間で合成して前記終点におけるバルク位相エラーを定める。

Description

信号処理装置、信号処理方法および信号処理プログラム
　本発明は、信号処理装置、信号処理方法および信号処理プログラムに関する。
　本願は、２０１９年１１月１４日に日本に出願された特願２０１９－２０６５０６号について優先権を主張し、それらの内容をここに援用する。
　光干渉断層撮影（ＯＣＴ：Optical Coherence Tomography）は、光の干渉性（コヒーレンス）を利用して試料（主に生体）の断層画像を取得する技術である。ＯＣＴによれば、試料の表面に限らず、その内部の構造を高い空間分解能で表す画像を取得することができる。ＯＣＴ技術の発達により、より高速かつ分解能の高い画像の撮影が可能となり、その機能の高度化が図られている。例えば、非特許文献１には、フーリエ領域ＯＣＴ（ＦＤ：Fourier Domain－ＯＣＴ）におけるディジタルデフォーカス（digital defocus）について記載されている。
　一般に、ＯＣＴによって得られた画像には、計測中の試料の動きや光が伝搬する大気や温度などの環境のゆらぎにより、バルク位相エラー（ＢＰＥ：Bulk Phase Error）が生じる。バルク位相エラーは、複数回の深さ方向への信号取得（Ａ－スキャン）の間で生ずる位相の揺らぎを指す。そのため、バルク位相エラーは、機能の高度化を阻む一因となりうる。非特許文献１に記載の手法は、バルク位相エラーが無視できるほど計測時間が短い一枚の断層画像（Ｂ－スキャン）を用いた近似的な手法に過ぎない。
　非特許文献２には、三次元のＯＣＴ画像（三次元ＯＣＴボリューム）に対して位相差法を用いて位相補正を行ってディジタルリフォーカスを行う手法について記載されている。位相差法は、三次元ＯＣＴボリュームの正面における２方向（（ｘ，ｙ）方向）の微分位相フィールドを計算し、微分位相フィールドを累積（cumulative summation）してバルク位相エラーを推定する手法である。
Yasuno et al., "Non-iterative numerical method for laterally superresolving Fourier domain optical coherence tomography,"OPTICS EXPRESS Vol.14, No. 3, p.1006-1020, February 6, 2006 Shemonski et al., "Three-dimensional motion correction using speckle and phase for in vivo computed optical interferometric tomography,"BIOMEDICAL OPTICS EXPRESS Vol.5, No.12, p.4131-4143, December 1, 2014.
　しかしながら、位相差法によれば、微分位相フィールドを積分する過程において誤差が累積するので、最終的に推定されるバルク位相エラーの精度が低下しがちであった。非特許文献２では、推定されたバルク位相エラーに対してスムージング処理を施して誤差を低減することを試みていた。この手法は、誤差の低減に対する理論的な根拠に乏しいため、むしろ推定精度の低下要因やバルク位相エラーの推定の失敗要因となることもあった。
　本発明は上記の点に鑑みてなされたものであり、本発明の課題は、バルク位相エラーの推定精度を向上することができる信号処理装置、信号処理方法および信号処理プログラムを提供することである。
（１）本発明は上記の課題を解決するためになされたものであり、本発明の一態様は、試料の状態を表す光干渉断層信号の光の照射方向に交差する平面における位相勾配を前記平面に配列されたサンプル点ごとに算出する位相勾配算出部と、バルク位相エラーが決定済のサンプル点である起点からバルク位相エラーが未決定のサンプル点である終点までの複数の経路のそれぞれについて、サンプル点ごとの前記位相勾配をそれぞれの経路に沿って積分して前記終点における経路別位相エラーを算出し、前記経路別位相エラーを前記複数の経路間で合成して前記終点におけるバルク位相エラーを定めるバルク位相エラー算出部と、を備える信号処理装置である。
（２）本発明の他の態様は、（１）の信号処理装置であって、前記バルク位相エラー算出部は、　前記経路別位相エラーの経路間の平均値を前記バルク位相エラーとして算出してもよい。
（３）本発明の他の態様は、（２）の信号処理装置であって、前記バルク位相エラー算出部は、前記平均値として、前記複数の経路それぞれの経路長に基づく重み係数を用いた前記経路別位相エラーの加重平均値を算出してもよい。
（４）本発明の他の態様は、（１）から（３）のいずれかの信号処理装置であって、前記複数の経路のそれぞれは、前記終点から所定の範囲内に含まれてもよい。
（５）本発明の他の態様は、（１）から（４）のいずれかの信号処理装置であって、前記複数の経路のそれぞれは、第ｎ（ｎは、１以上かつＮ－１以下の整数）サンプル点から第ｎ＋１サンプル点の区間である第ｎ素片（Ｎは、２以上の整数）を含み、前記第ｎ＋１サンプル点は、前記第ｎサンプル点から前記光干渉断層信号を取得した光学系の空間分解能以下の範囲内のサンプル点であってもよい。
（６）本発明の他の態様は、（１）から（５）のいずれかの信号処理装置であって、前記起点と前記終点の位置関係に対応付けて前記複数の経路を示す経路パターンを予め複数通り設定しておき、バルク位相エラーの算出対象とする対象サンプル点を終点として、前記位置関係を満たす経路パターンで示される複数の経路を選択してもよい。
（７）本発明の他の態様は、（６）の信号処理装置であって、前記バルク位相エラー算出部には、前記平面における１つのサンプル点を基準点とし、前記基準点に近接するサンプル点ほど優先されるように、バルク位相エラーを算出する順序を予め定めておき、前記順序が第ｍ回（ｍは、１以上の整数）のサンプル点を前記対象サンプル点とするとき、前記順序が第ｍ－１回のサンプル点を起点として、前記位置関係を満たす経路パターンで示される前記複数の経路を選択してもよい。
（８）本発明の他の態様は、（１）から（７）のいずれかの信号処理装置であって、前記終点におけるバルク位相エラーを、前記終点と前記平面に射影された位置が共通する三次元空間内に配列されたサンプル点における前記光干渉断層信号の信号値の位相から減算する信号補正部を備えてもよい。
（９）本発明の他の態様は、信号処理装置における方法であって、試料の状態を表す光干渉断層信号の光の照射方向に交差する平面における位相勾配を前記平面に配列されたサンプル点ごとに算出する位相勾配算出過程と、バルク位相エラーが決定済のサンプル点である起点からバルク位相エラーが未決定のサンプル点である終点までの複数の経路のそれぞれについて、サンプル点ごとの前記位相勾配をそれぞれの経路に沿って積分して前記終点における経路別位相エラーを算出し、前記経路別位相エラーを前記複数の経路間で合成して前記終点におけるバルク位相エラーを定めるバルク位相エラー算出過程と、を有する信号処理方法である。
（１０）本発明の他の態様は、試料の状態を表す光干渉断層信号の光の照射方向に交差する平面における位相勾配を前記平面に配列されたサンプル点ごとに算出する位相勾配算出手順と、バルク位相エラーが決定済のサンプル点である起点からバルク位相エラーが未決定のサンプル点である終点までの複数の経路のそれぞれについて、サンプル点ごとの前記位相勾配をそれぞれの経路に沿って積分して前記終点における経路別位相エラーを算出し、前記経路別位相エラーを前記複数の経路間で合成して前記終点におけるバルク位相エラーを定めるバルク位相エラー手順と、を実行させるための信号処理プログラムである。
　本発明によれば、バルク位相エラーの推定精度を向上することができる。
本実施形態に係る光干渉断層計の一例を示す構成図である。 本実施形態に係る信号処理装置の機能構成例を示すブロック図である。 本実施形態に係るバルク位相エラーの演算処理の概要を示す説明図である。 本実施形態に係る積分経路の例を示す説明図である。 本実施形態に係る経路パターンの第1例を示す図である。 本実施形態に係る経路パターンの第２例を示す図である。 本実施形態に係る経路パターンの第３例を示す図である。 本実施形態に係る経路パターンの第４例を示す図である。 本実施形態に係る経路パターンの第５例を示す図である。 本実施形態に係る経路パターンの第６例を示す図である。 本実施形態に係る経路パターンの第７例を示す図である。 本実施形態に係る経路パターンの第８例を示す図である。 本実施形態に係る経路パターンの第９例を示す図である。 本実施形態に係る経路パターンの起点と終点のセットの第1例を示す図である。 本実施形態に係る経路パターンの起点と終点のセットの第２例を示す図である。 本実施形態に係る経路パターンの起点と終点のセットの第３例を示す図である。 本実施形態に係る経路パターンの起点と終点のセットの第４例を示す図である。 本実施形態に係る経路パターンの起点と終点のセットの第５例を示す図である。 本実施形態に係る経路パターンの起点と終点のセットの第６例を示す図である。 本実施形態に係る経路パターンの起点と終点のセットの第７例を示す図である。 本実施形態に係るバルク位相エラーの繰り返し演算の初期状態の例を示す図である。 本実施形態に係るバルク位相エラーの繰り返し演算の１回目の例を示す図である。 本実施形態に係るバルク位相エラーの繰り返し演算の２回目の例を示す図である。 本実施形態に係るバルク位相エラーの繰り返し演算の３回目の例を示す図である。 本実施形態に係るバルク位相エラーの繰り返し演算の他の例を示す図である。 ＯＣＴ画像の例を示す図である。 位相差法に基づくバルク位相エラーの例を示す図である。 本実施形態に係るバルク位相エラーの例を示す図である。 ＯＣＴ画像のスペクトルの例を示す図である。 位相差法に基づくバルク位相エラーを用いて補正されたＯＣＴ画像のスペクトルの例を示す図である。 本実施形態に係るバルク位相エラーを用いて補正されたＯＣＴ画像のスペクトルの例を示す図である。 空間スペクトルのｘ方向平均値の例を示す図である。 空間スペクトルのｙ方向平均値の例を示す図である。 空間スペクトルのｘ方向標準偏差の例を示す図である。 空間スペクトルのｙ方向標準偏差の例を示す図である。
　以下、図面を参照しながら本発明の実施形態について説明する。
　図１は、本実施形態に係る光干渉断層計１の一例を示す構成図である。
　光干渉断層計１は、ＯＣＴを用いて試料の状態を観測するための観測システムを構成する。
　光干渉断層計１は、試料Ｓｍに光を照射し、試料Ｓｍから反射した反射光と、参照鏡４０（後述）で反射された参照光とを干渉させて生じた干渉光を取得し、取得した干渉光から試料Ｓｍの表面とその内部の状態を示す画像を生成する装置である。
　試料Ｓｍとする観測対象の物体は、例えば、人間もしくは動物の生体、非生物のいずれでもよい。生体は、眼底、血管、歯牙、皮下組織などでありうる。非生物は、電子部品、機械部品など人工的な構造体、石材、鉱物などの天然の構造体、特定の形状を有しない物質のいずれでもよい。
　光干渉断層計１は、光源１０、ビームスプリッタ２０、コリメータ３０ａ、３０ｂ、５０ａ、５０ｂ、参照鏡４０、ガルバノミラー６０ａ、６０ｂ、分光器７０、および信号処理装置１００を含んで構成される。これらの構成部材のうち、ビームスプリッタ２０、コリメータ３０ａ、３０ｂ、５０ａ、５０ｂ、参照鏡４０、ガルバノミラー６０ａ、６０ｂ、および分光器７０は、干渉計と呼ばれる光学系を構成する。図１に例示される干渉計は、光ファイバＦを備えるマイケルソン干渉計である。より具体的には、光源１０と分光器７０、コリメータ３０ａとコリメータ５０ａは、それぞれ光ファイバＦを用いてビームスプリッタに接続されている。光ファイバＦは、光源１０から照射される光の波長帯を含む伝送帯域を有する。
　光干渉断層計１は、例えば、フーリエドメインＯＣＴ（ＦＤ－ＯＣＴ：Fourier-Domain　ＯＣＴ）である。ＦＤ－ＯＣＴは、例えば、スペクトラルドメインＯＣＴ（ＳＤ－ＯＣＴ：Spectral-DomainＯＣＴ）、波長掃引ＯＣＴ（Swept-Source ＯＣＴ）等のいずれであってもよいが、これには限られない。光干渉断層計１は、例えば、タイムドメインＯＣＴ（Time-Domain ＯＣＴ：ＴＤ－ＯＣＴ）であってもよい。光干渉断層計１がＴＤ－ＯＣＴである場合には、参照鏡４０の位置は固定されず、光源１０から参照鏡４０までの光路長を可変とする駆動機構を備える。
　光源１０は、例えば、近赤外の波長（例えば、８００～１０００ｎｍ）を有するプローブ光を照射する。光源１０は、例えば、超短波パルスレーザ、ＳＬＤ（スーパールミネセントダイオード；Superluminescent Diode）などの波長掃引光源である。光源１０から照射された光は、光ファイバＦ内部で導光され、ビームスプリッタ２０に入射する。
　ビームスプリッタ２０は、入射した光をコリメータ３０ａに向けて導光される光（以下、参照光）と、コリメータ５０ａに向けて導光される光（以下、測定光）に分離する。ビームスプリッタ２０は、例えば、キューブビームスプリッタなどである。
　コリメータ３０ａは、ビームスプリッタ２０から導光される参照光を平行光に変化させ、平行光をコリメータ３０ｂに向けて出射する。
　コリメータ３０ｂは、コリメータ３０ａから入射される平行光を集光し、集光した参照光を参照鏡４０に向けて出射する。なお、コリメータ３０ｂは、参照鏡４０で反射した参照光を入射し、平行光に変化させ、平行光をコリメータ３０ａに向けて出射する。
　コリメータ３０ａは、コリメータ３０ｂから入射した平行光を集光し、ビームスプリッタ２０に向けて導光する。
　他方、コリメータ５０ａは、ビームスプリッタ２０から導光される測定光を平行光に変化させ、平行光をガルバノミラー６０ａに向けて出射する。ガルバノミラー６０ａ、６０ｂの表面において、コリメータ５０ａから入射される平行光は、それぞれ反射され、コリメータ５０ｂに向けて出射される。コリメータ５０ｂは、コリメータ５０ａからガルバノミラー６０ａ、６０ｂを経由して入射された平行光を集光し、集光した測定光を試料Ｓｍに照射する。試料Ｓｍに照射される測定光は、試料Ｓｍの反射面において反射されコリメータ５０ｂに入射される。反射面は、例えば、試料Ｓｍと試料Ｓｍの周囲環境（例えば、大気）との境界面に限られず、試料Ｓｍ内部における屈折率が異なる材質間もしくは組織間の境界面となりうる。以下、試料Ｓｍの反射面において反射され、コリメータ５０ｂに入射される光を反射光と呼ぶ。
　コリメータ５０ｂは、入射された反射光をガルバノミラー６０ｂに向けて出射する。ガルバノミラー６０ｂ、６０ａそれぞれの表面において、それぞれ反射され、コリメータ５０ａに向けて出射される。コリメータ５０ａは、コリメータ５０ａからガルバノミラー６０ａ、６０ｂを経由して入射された平行光を集光し、集光した反射光をビームスプリッタ２０に向けて導光する。
　ビームスプリッタ２０は、参照鏡４０で反射された参照光と試料Ｓｍで反射された反射光とを、光ファイバＦを経由して分光器７０に導光する。
　分光器７０は、その内部に回折格子と受光素子を備える。回折格子は、ビームスプリッタ２０から導光された参照光と反射光を分光する。分光した参照光と反射光は、互いに干渉し、干渉した干渉光となる。受光素子は、干渉光を検出し、検出した干渉光に基づく信号（以下、検出信号）を生成する。受光素子は、所定の撮像面に複数の画素を二次元配置して構成される。受光素子は、生成した検出信号を信号処理装置１００に出力する。
　信号処理装置１００は、分光器７０から入力された検出信号に基づいて試料Ｓｍの状態を表すＯＣＴ信号を取得する。信号処理装置１００は、観測点を所定の順序に従って順次変更し、個々の観測点において、または所定の単位をなす複数の観測点ごとに一括して取得された検出信号を、観測点間で順次累積して所定の領域内のＯＣＴ信号を生成する。信号処理装置１００は、サンプル点ごとの補正後の信号値を画素ごとの輝度値に変換し、個々のサンプル点にそれぞれ対応する画素ごとの輝度値を示す画像データを生成する。
　本実施形態において、信号処理装置１００は、取得したＯＣＴ信号に対して画像データを生成する前に次の処理を行う。信号処理装置１００は、所定の平面に配列されたサンプル点ごとにＯＣＴ信号の位相勾配を算出する。信号処理装置１００は、バルク位相エラーを決定済みのサンプル点を起点とし、バルク位相エラーが未決定のサンプル点を終点と定め、起点からある１つの終点までの複数の経路を設定する。そして信号処理装置１００は、設定した経路ごとに、その経路上のサンプル点ごとの位相勾配を起点から終点までそれぞれの経路に沿って積分してバルク位相エラー（以下、経路別位相エラー）を算出する。信号処理装置１００は、経路別位相誤差を設定した経路間で合成し、終点におけるバルク位相エラーを定める。信号処理装置１００は、その平面における終点とするサンプル点のそれぞれについてバルク位相エラーを定める処理を実行する。バルク位相エラーとは、試料の動きや光が伝搬する大気や温度のゆらぎを原因としてＯＣＴ信号の位相に含まれるエラーを意味する。
　信号処理装置１００は、ＯＣＴ信号を示すサンプル点ごとの信号値の位相からバルク位相エラーを差し引いて、その信号値の位相を補正する。
　次に、本実施形態に係る信号処理装置１００の機能構成例について説明する。
　図２は、本実施形態に係る信号処理装置１００の機能構成例を示すブロック図である。
　信号処理装置１００は、制御部１１０と記憶部１８０を含んで構成される。制御部１１０は、例えば、ＣＰＵ（Central Processing Unit）等のプロセッサが予め記憶部１８０に記憶させておいたプログラムを読み出し、読み出したプログラムに記述された指令で指示される処理を行って、その機能を奏する。本願では、プログラムに記述された指令で指示される処理を行うことを、プログラムを実行する、プログラムの実行、などと呼ぶことがある。制御部１１０の一部または全部は、ＬＳＩ（Large Scale Integration）、ＡＳＩＣ（Application Specific Integrated Circuit）等のハードウェアもしくは専用のハードウェアで構成されてもよい。
　制御部１１０は、光学系制御部１２０、検出信号取得部１３０、位相勾配算出部１４０、バルク位相エラー算出部１５０、信号補正部１６０、および画像処理部１７０を含んで構成される。
　光学系制御部１２０は、ガルバノミラー６０ａ、６０ｂの位置を可変にする駆動機構を駆動させ、試料Ｓｍの観測点であるサンプル点を走査する。試料Ｓｍのサンプル点は、試料Ｓｍの深さ方向と、その方向に交差する方向に交差する方向（例えば、垂直方向）のそれぞれについて走査される。以下の説明では、試料の深さ方向を、三次元直交座標系におけるｚ方向と呼び、ｚ方向と相互に直交する第１方向、第２方向をそれぞれｘ方向、ｙ方向と呼ぶことがある。試料の深さ方向への信号取得は、Ａ－スキャン（A-scan）と呼ばれ、その方向に交差する方向（例えば、ｘ方向またはｙ方向）への信号取得は、Ｂ－スキャン（B-scan）と呼ばれることがある。
　検出信号取得部１３０は、分光器７０から検出信号を逐次に取得する。検出信号取得部１３０は、取得した検出信号に基づいて試料Ｓｍの深さ方向の測定光の反射された光（反射光）の分布を信号値として試料Ｓｍが存在する三次元空間において所定の間隔で配列されたサンプル点ごとに算出する。試料Ｓｍの深さ方向は、測定光の入射方向（ｚ方向）に相当する。算出される屈折率分布は、反射光の強度分布に相当する。検出信号取得部１３０は、ｚ方向に交差する面（例えば、ｘ－ｙ平面）に沿った信号取得により変更される観測点ごとに反射光の強度分布を算出する処理を繰り返す。これにより、検出信号取得部１３０は、観測可能とする三次元領域（以下、観測可能領域）内の反射光の分布を表すデータ（以下、ＯＣＴ信号）を取得することができる。ＯＣＴ信号は、観測可能領域において一定の間隔で分布したサンプル点ごとの信号値を示す。ＯＣＴ信号は、観測可能領域における試料Ｓｍの状態を表す三次元のデータであり、ＯＣＴボリューム（OCT volume）、三次元ＯＣＴ信号（図３）などとも呼ばれる。検出信号取得部１３０は、取得したＯＣＴ信号を記憶部１８０に記憶する。ｚ方向に交差する面は、必ずしもｚ方向に直交していなくてもよく、ｚ方向と平行でなければよい。また、ｚ方向に交差する二次元の面において、個々のサンプル点は、必ずしも、直交格子の各格子点上に配置されていなくてもよい。個々のサンプル点は、例えば、斜め格子の各格子点上に配置されてもよいし、非周期的に配置されてもよい。なお、隣接するサンプル点の間隔は、光学系の空間分解能以下とする。言い換えれば、隣接するサンプル点間において試料の構造に起因する位相の成分（試料位相：sample phase）に有意差が生じない。後述するように、位相勾配の要素として隣接するサンプル点間での信号値の位相差を算出することにより、試料位相の成分が相殺され、バルク位相エラーの成分が残される。
　上記の説明では、対象とするサンプル点の変更を伴う信号取得もしくは複数のサンプル点に対する一括した信号取得に対して「スキャン」と呼称しているが、必ずしも光学系を構成する部材（例えば、コリメータ５０ｂ、試料Ｓｍの支持台）または検出器の駆動による走査を伴わない場合も含まれうる。その場合には、信号取得に係る機能を検出信号取得部１３０が担い、光学系制御部１２０において省略される。例えば、ＯＣＴの方式としてＦＤ－ＯＣＴを用いてｘ－ｙ平面上のある１点においてＡ－スキャンを行う際、検出信号取得部１３０は、試料Ｓｍの深さにより異なる周波数成分を有する反射光に基づく干渉光を検出信号として検出する。検出信号取得部１３０は、検出信号をフーリエ変換して周波数ごとの変換係数を算出し、周波数に対応付けられた深さの変換係数をその深さにおける信号値として定めることができる。さらに、ｘ方向またはｙ方向への駆動による走査を要しない場合には、光学系制御部１２０が省略されうる。
　位相勾配算出部１４０は、記憶部１８０からＯＣＴ信号を読み取り、読み取ったＯＣＴ信号が示すサンプル点ごとの信号値に対して公知の手法を用いて位相アンラップ（phase unwrapping）を行い、その信号値の位相を調整して調整済の位相（unwrapped phase）を算出する。位相アンラップは、－πからπの値域に制限された位相について空間的な連続性を担保するように、その制限を解除する処理である。位相勾配算出部１４０は、サンプル点ごとの調整済の位相に基づいて、所定の二次元平面上に配列されるサンプル点ごとに位相勾配を算出する（図３：位相勾配解析S01）。位相勾配は、その二次元平面に平行で互いに直交する第１の方向（例えば、ｘ方向）と第２の方向（例えば、ｙ方向）のそれぞれに対する偏微分を要素として含む二次元のベクトルである。二次元平面は、ｚ方向に交差する面であれば、ｚ方向に直交する面であっても、直交しない面であってもよい。また、二次元平面は、三次元の観測可能領域内の任意の点から射影可能とする平面であれば観測可能領域を横断する平面であってもよいし、仮想的な平面であってもよい。二次元平面は、例えば、三次元の観測可能領域内に配列されたサンプル点は、その二次元平面に垂直な第３の方向（例えば、ｚ方向）に射影した点と対応付けられる。ｚ方向に交差する二次元平面内のサンプル点ごとに位相勾配を算出するのは、バルク位相エラーがｚ方向には依存せず、主にその二次元平面内の座標に依存するためである。なお、以下の説明では二次元平面上に射影されたサンプル点が分布している領域を二次元領域と呼ぶことがある。
　位相勾配算出部１４０は、二次元領域におけるサンプル点ごとに算出した位相勾配を示すデータを位相勾配データ（図３：二次元位相勾配Pg01）として記憶部１８０に記憶する。
　バルク位相エラー算出部１５０は、記憶部１８０から位相勾配データを読み取り、読み取った位相勾配データが示すサンプル点ごとの位相勾配に基づいて、以下に説明する処理を行って、二次元領域に配列されたそれぞれのサンプル点におけるバルク位相エラーを算出する。
　バルク位相エラー算出部１５０は、バルク位相エラーを決定済のサンプル点の１つまたは複数個のそれぞれを起点とし、起点からバルク位相エラーが未決定のサンプル点の１つである終点までの複数の経路を設定する。バルク位相エラー算出部１５０は、設定した経路のそれぞれについて、サンプル点ごとの位相勾配を、それぞれの経路に沿って起点から終点まで積分（図３：経路積分による位相推定S02）して終点における積分値を経路別位相エラーとして算出する。そして、バルク位相エラー算出部１５０は、経路ごとに算出した経路別位相エラーを設定した複数の経路間で合成して算出される合成値を終点におけるバルク位相エラー（図３：バルク位相エラーBp01）として定める。但し、バルク位相エラー算出部１５０は、二次元平面上のサンプル点のうちの１点を基準点（例えば、二次元領域の原点）とし、基準点におけるバルク位相エラーを予め所定の基準値（例えば、０）として予め定めておく。バルク位相エラー算出部１５０は、定めたサンプル点ごとのバルク位相エラーを示すバルク位相エラーデータを記憶部１８０に記憶する。
　信号補正部１６０は、記憶部１８０からバルク位相エラーデータとＯＣＴ信号を読み出す。信号補正部１６０は、読み出したバルク位相エラーデータが示すサンプル点ごとのバルク位相エラーに基づいて、ＯＣＴ信号が示すサンプル点ごとの信号値の位相を補正する処理を行う。
　信号補正部１６０は、バルク位相エラーデータが示す二次元領域上のサンプル点（以下、二次元サンプル点）の位置と、二次元平面に射影された点の位置が共通する観測可能領域内のサンプル点（以下、三次元サンプル点）を対応する三次元サンプル点として特定する。特定される三次元サンプル点の数は、１つの二次元サンプル点に対して複数となりうる。信号補正部１６０は、特定した三次元サンプル点における信号値の位相を、対応する二次元サンプル点のバルク位相エラーを減算することにより補正する。信号補正部１６０は、三次元サンプル点ごとに位相を補正した信号値に対して、さらに公知の補正処理（例えば、再焦点化（リフォーカシング（refocusing））、スペックル抑圧など）を行ってもよい。信号補正部１６０は、補正処理により得られた三次元サンプル点ごとの信号値を示す補正ＯＣＴ信号を記憶部１８０に記憶する。
　画像処理部１７０は、記憶部１８０から補正ＯＣＴ信号を読み出し、三次元サンプル点ごとの信号値を所定の変換関数を用いて輝度値に変換する。変換される輝度値は、画素ごとのビット深度で表現可能な値域内の値をとる。観察対象の二次元平面（以下、観察対象平面）におけるサンプル点ごとの輝度値を画素ごとの輝度値として有する画像データを生成する。本願では、サンプル点を画素と呼ぶことがある。
　画像処理部１７０は、操作入力部（図示せず）から入力された操作信号に基づいて観測可能領域の少なくとも一部を横切る断面を観察対象の二次元平面として定めてもよい。操作信号により、観察対象平面の要件として、例えば、試料表面からの深度、試料表面に対する観察方向などが指示される。
　操作入力部は、例えば、ボタン、つまみ、ダイヤル、マウス、ジョイスティックなど、ユーザの操作を受け付け、受け付けた操作に応じた操作信号を生成する部材を含んで構成されてもよい。操作入力部は、他の機器（例えば、リモートコントローラ等の携帯機器）から操作信号を無線または有線で受信する入力インタフェースであってもよい。
　画像処理部１７０は、生成した画像データを表示部（図示せず）に出力し、画像データに基づくＯＣＴ画像を表示させてもよい。表示部は、例えば、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、有機エレクトロルミネッセンスディスプレイ（ＯＬＥＤ：Organic Electroluminescence Display）などのいずれであってもよい。
　画像処理部１７０は、生成した画像データを記憶部１８０に記憶してもよいし、通信部（図示せず）を経由して、他の機器に無線または有線で送信してもよい。
　記憶部１８０は、上記のプログラムの他、制御部１１０が実行する処理に用いられる各種のデータ、制御部１１０が取得した各種のデータを記憶する。
　記憶部１８０は、例えば、ＲＯＭ（Read Only Memory）、フラッシュメモリ、ＨＤＤ（Hard Disk Drive）などの不揮発性の（非一時的）記憶媒体を含んで構成される。記憶部１８０は、例えば、ＲＡＭ（Random Access Memory）、レジスタなどの揮発性の記憶媒体を含んで構成される。
（位相勾配）
　次に、位相勾配の算出例について説明する。位相勾配算出部１４０は、三次元サンプル点（ｘ，ｙ，ｚ）ごとの調整済の位相φ（ｘ，ｙ，ｚ）に基づいて、二次元領域に配列される二次元サンプル点（ｘ，ｙ）ごとに位相勾配（Ｂｘ，Ｂｙ）を算出する。図３に示す例では、三次元サンプル点（ｘ，ｙ，ｚ）は直方体の観測可能領域内に、ｘ方向、ｙ方向、ｚ方向にそれぞれ所定の間隔で分布している。ｘ方向、ｙ方向、ｚ方向は、互いに直交し、観測可能領域の外縁をなす各辺の方向に相当する。二次元領域は、ｘ方向とｙ方向に平行な二次元のｘ－ｙ平面である。三次元サンプル点（ｘ，ｙ，ｚ）を二次元領域上に射影した二次元サンプル点（ｘ，ｙ）の座標は、三次元サンプル点（ｘ，ｙ，ｚ）のうち、ｚ方向の座標値を無視して得ることができる。言い換えれば、二次元サンプル点（ｘ，ｙ）に対応する三次元サンプル点（ｘ，ｙ，ｚ）は、その座標が、二次元サンプル点（ｘ，ｙ）のｘ座標値、ｙ座標値とそれぞれ共通なｘ座標値、ｙ座標値で示されるサンプル点である。
　位相勾配算出部１４０は、例えば、式（１）に示すように三次元サンプル点（ｘ，ｙ，ｚ）における信号値のスペクトルＡ（ｘ，ｙ，ｚ）とｘ方向に隣接する三次元サンプル点（ｘ＋Δｘ，ｙ，ｚ）における信号値のスペクトルの複素共役Ａ＊（ｘ＋Δｘ，ｙ，ｚ）との乗算値の観測領域内のｚ方向にわたる総和の偏角を、対応する二次元サンプル点（ｘ，ｙ）における位相勾配のｘ方向成分Ｂｘ（ｘ，ｙ）として算出する。
Figure JPOXMLDOC01-appb-M000001
　位相勾配算出部１４０は、対応する二次元サンプル点（ｘ，ｙ）における位相勾配のｙ方向成分Ｂｙ（ｘ，ｙ）についても、式（２）に示すように三次元サンプル点（ｘ，ｙ，ｚ）における信号値のスペクトルＡ（ｘ，ｙ，ｚ）とｙ方向に隣接する三次元サンプル点（ｘ，ｙ＋Δｙ，ｚ）における信号値のスペクトルの複素共役Ａ＊（ｘ，ｙ＋Δｙ，ｚ）に基づいて算出することができる。
Figure JPOXMLDOC01-appb-M000002
　式（１）、（２）において、Δｘ、Δｙは、それぞれｘ、ｙ方向への三次元サンプル点の間隔を示すとともに、二次元サンプル点の間隔を示す。ａｒｇ（…）は、複素数…の偏角を示す。
　式（１）、（２）に示すように、ｘ－ｙ平面に直交するｚ方向に分布するサンプル点間の乗算値に対する偏角を求めることで、個々の二次元サンプル点間に生ずるランダムな試料位相の残差が累積されない反面、バルク位相エラーが累積される。そのため、試料位相の残差に起因する誤差を低減することができる。また、強度が高いサンプル点が、ノイズの影響が顕著となりがちな強度が低いサンプル点よりも重視して位相勾配が算出されるため、ノイズの影響が低減される。二次元サンプル点ごとに算出された位相勾配Ｂは、上記のようにバルク位相エラーの算出に用いられる。
　なお、後述するように積分経路は、必ずしもｘ方向またはｙ方向に平行な部分だけに限らず、ｘ方向とｙ方向の両者に交差する方向をとる部分を含みうる。その場合には、位相勾配算出部１４０は、式（３）に基づいて、その部分のサンプル点（ｘ，ｙ）における、次のサンプル点（ｘ＋Δｘ、ｙ＋Δｙ）と位相勾配Ｂｘｙ（ｘ，ｙ）を算出することができる。
Figure JPOXMLDOC01-appb-M000003
（経路積分）
　次に、位相勾配の経路積分の例について説明する。
　バルク位相エラー算出部１５０には、位相エラーの決定対象とする二次元サンプル点（以下、位相決定対象画素とも呼ぶ）ごとに、個々の位相決定対象画素を終点とする複数の積分経路を予め設定しておく。個々の積分経路の起点は、位相エラーを決定済の二次元サンプル点（以下、位相決定済画素とも呼ぶ）とする。但し、起点は、終点から所定の範囲内に所在する二次元サンプル点とする。これにより、積分経路の数が制限され過大にならない。
　バルク位相エラー算出部１５０は、例えば、式（４）に示すように経路Ｐｔ上のサンプル点（ｘ，ｙ）ごとの位相勾配Ｂを、それぞれの経路Ｐｔの起点Ｏｐから終点Ｄｐまで経路積分を行って得られる終点における積分値を経路別位相エラーφｅ，Ｐｔとして算出する。
Figure JPOXMLDOC01-appb-M000004
　式（４）において、Ｂ（ｘ，ｙ）は、位相勾配を示し、サンプル点（ｘ，ｙ）におけるｘ方向成分Ｂｘ（ｘ，ｙ）とｙ方向成分Ｂｙ（ｘ，ｙ）を要素として含む二次元のベクトルである。ｎ（ｘ，ｙ）は、経路Ｐｔの接線方向の方向ベクトルを示す。より具体的には、方向ベクトルｎ（ｘ，ｙ）として、サンプル点（ｘ，ｙ）を起点とし、経路Ｐｔ上に隣接するサンプル点を終点とするベクトルが利用可能である。・は、ベクトルの内積を示す。
　そして、バルク位相エラー算出部１５０は、経路Ｐｔについて算出された経路別位相エラーφｅ，Ｐｔを、終点Ｄｐが共通な経路間で合成し、その終点Ｄｐにおけるバルク位相エラーφｅを算出する。
　次に、二次元平面上に設定される積分経路の例について説明する。
　図４は、本実施形態に係る積分経路の例を示す説明図である。図４に示す例では、起点とする位相決定済画素Ｏｐ００と終点とする位相決定対象画素Ｄｐ１１に対して、４つの経路Ｐｔ０１－Ｐｔ０４が設定されている。位相決定済画素Ｏｐ００、位相決定対象画素Ｄｐ１１の二次元平面上の座標は、それぞれ（０，０）、（１，１）である。終点は、起点に対してｘ座標とｙ座標の双方に交差する斜め方向に隣接している。ｘ座標、ｙ座標が、それぞれ間隔Δｘ、Δｙが１となるように正規化されている。経路Ｐｔ０１－Ｐｔ０４は、それぞれ１つのサンプル点を起点とし、ｘ方向またはｙ方向に隣接するサンプル点を終点とする素片を順次接続して形成される。個々の素片を示すベクトルが式（４）に示す方向ベクトルとして経路別位相エラーの算出に用いられる。この例では、個々の素片の長さ（以下、素片長）が１であるため、バルク位相エラー算出部１５０は、それぞれの経路Ｐｔ０１－Ｐｔ０４を形成する素片の数を経路の長さ（以下、経路長）として定めることができる。経路Ｐｔ０１、Ｐｔ０２、Ｐｔ０３、Ｐｔ０４の経路長は、それぞれ２、２、４、４となる。
　従って、バルク位相エラー算出部１５０は、位相決定対象画素Ｄｐ１１における経路Ｐｔ０１、Ｐｔ０２、Ｐｔ０３、Ｐｔ０４に係る経路別位相エラーφｅ，Ｐｔ０１、φｅ，Ｐｔ０２、φｅ，Ｐｔ０３、φｅ，Ｐｔ０４を、それぞれＢｘ（０，０）＋Ｂｙ（１，０）、Ｂｙ（０，０）＋Ｂｘ（０，１）、Ｂｘ（０，０）＋Ｂｘ（１，０）＋Ｂｙ（２，０）－Ｂｘ（２，１）、Ｂｙ（０，０）＋Ｂｙ（０，１）＋Ｂｘ（０，２）－Ｂｙ（１，２）と算出することができる。
　そして、バルク位相エラー算出部１５０は、経路別位相エラーを経路間で合成したバルク位相エラーとして、例えば、経路間の単純平均値を算出する。図４に示す例では、バルク位相エラーを算出する過程で、経路別位相エラーの総和を経路数である４を除算して正規化する処理が含まれる。
　但し、経路長に応じて位相勾配に重畳される誤差の影響が異なりうる。そこで、バルク位相エラー算出部１５０は、経路ごとの重み係数として経路長に依存する値を用いて、経路別位相エラーの加重平均値をバルク位相エラーとして算出してもよい。
　例えば、バルク位相エラー算出部１５０は、経路ごとの重み係数として、ａｐ（ａは、０より大きく１より小さい実数、ｐは、経路長）に比例し、経路間の総和が１に正規化された実数を用いることができる。その場合、図４に示す例では、経路別位相エラーφｅ，Ｐｔ０１、φｅ，Ｐｔ０２、φｅ，Ｐｔ０３、φｅ，Ｐｔ０４に対する重み係数は、それぞれａ２／２（ａ２＋ａ４）、ａ２／２（ａ２＋ａ４）、ａ４／２（ａ２＋ａ４）、ａ４／２（ａ２＋ａ４）となる。
　バルク位相エラー算出部１５０は、経路ごとの重み係数として、経路長ｐに反比例し、経路間の総和が１に正規化された値を用いてもよい。その場合、図４に示す例では、経路別位相エラーφｅ，Ｐｔ０１、φｅ，Ｐｔ０２、φｅ，Ｐｔ０３、φｅ，Ｐｔ０４に対する重み係数は、それぞれ１／３、１／３、１／６、１／６となる。これらの手法によれば、経路長が大きい経路の経路別位相エラーのバルク位相エラーへの寄与を小さくすることができる。
　また、バルク位相エラー算出部１５０は、式（５）を用いて終点（ｘ０，ｙ０）におけるバルク位相エラーφｂ （ｔ）（ｘ０，ｙ０）を算出してもよい。式（５）において、φｅ，Ｐｔ（ｘ０，ｙ０）は、経路Ｐｔに対する経路別位相エラーを示す。式（５）は、大きさが１であり、経路別位相エラーを偏角とする単位ベクトルの経路間の総和である総和ベクトルを算出し、算出した総和ベクトルの偏角をバルク位相エラーφｂ （ｔ）（ｘ０，ｙ０）として定めることを示す。
Figure JPOXMLDOC01-appb-M000005
　なお、位相決定済画素と位相決定対象画素との位置関係に応じて、種々の積分経路が設定可能である。以下の説明では、位相決定済画素と位相決定対象画素との位置関係に対応する複数の積分経路の組を経路パターンと呼ぶ。図４に例示される経路パターンを経路パターンＰｎ００と呼ぶ。
　図５Ａ－図５Ｉは、本実施形態に係る経路パターンのその他の例を示す図である。但し、図５Ａ－図５Ｉのそれぞれに示す経路パターンＰｎ０１－Ｐｎ０９は、いずれも終点が位相決定対象画素Ｄｐ１１であり、起点とする位相決定済画素の位置は、位相決定対象画素Ｄｐ１１に対してｘ軸の負方向とｙ軸の負方向の一方または両方に隣接する。起点とする位相決定済画素の数は、１個に限られず、２個以上となりうる。経路パターンＰｎ００－Ｐｎ０２では、起点とする位相決定済画素の数は１個である。これに対し、経路パターンＰｎ０３－Ｐｎ０５では、起点とする位相決定済画素の数は２個である。経路パターンＰｎ０６－Ｐｎ０９では、起点とする位相決定済画素の数は３個である。
　図５Ａに例示する経路パターンＰｎ０１は、３つの経路Ｐｔ１１－Ｐｔ１３を有する。経路Ｐｔ１１－Ｐｔ１３は、それぞれ位相決定済画素Ｏｐ０１を起点とする。経路Ｐｔ１１、Ｐｔ１２、Ｐｔ１３の経路長は、それぞれ３、１、３となる。
　図５Ｂに例示する経路パターンＰｎ０２は、３つの経路Ｐｔ２１－Ｐｔ２３を有する。経路Ｐｔ２１－Ｐｔ２３は、それぞれ位相決定済画素Ｏｐ１０を起点とする。経路Ｐｔ２１、Ｐｔ２２、Ｐｔ２３の経路長は、それぞれ３、１、３となる。
　図５Ｃに例示する経路パターンＰｎ０３は、２つの経路Ｐｔ３１、Ｐｔ３２を有する。経路Ｐｔ３１、Ｐｔ３２は、それぞれ位相決定済画素Ｏｐ０１、Ｏｐ０２を起点とする。経路Ｐｔ３１、Ｐｔ３２の経路長は、それぞれ１、３となる。
　図５Ｄに例示する経路パターンＰｎ０４は、２つの経路Ｐｔ４１、Ｐｔ４２を有する。経路Ｐｔ４１、Ｐｔ４２は、それぞれ位相決定済画素Ｏｐ１０、Ｏｐ２０を起点とする。経路Ｐｔ４１、Ｐｔ４２の経路長は、それぞれ１、３となる。
　図５Ｅに例示する経路パターンＰｎ０５は、４つの経路Ｐｔ５１－Ｐｔ５４を有する。経路Ｐｔ５１、Ｐｔ５２、Ｐｔ５３、Ｐｔ５４は、それぞれ位相決定済画素Ｏｐ０１、Ｏｐ１０、Ｏｐ０１、Ｏｐ１０を起点とする。経路Ｐｔ５１、Ｐｔ５２、Ｐｔ５３、Ｐｔ５４の経路長は、それぞれ１、１、３、３となる。
　図５Ｆに例示する経路パターンＰｎ０６は、３つの経路Ｐｔ６１－Ｐｔ６３を有する。経路Ｐｔ６１、Ｐｔ６２、Ｐｔ６３は、それぞれ位相決定済画素Ｏｐ００、Ｏｐ０１、Ｏｐ０２を起点とする。経路Ｐｔ６１、Ｐｔ６２、Ｐｔ６３の経路長は、それぞれ２、１、２となる。
　図５Ｇに例示する経路パターンＰｎ０７は、３つの経路Ｐｔ７１－Ｐｔ７３を有する。経路Ｐｔ７１、Ｐｔ７２、Ｐｔ７３は、それぞれ位相決定済画素Ｏｐ００、Ｏｐ１０、Ｏｐ２０を起点とする。経路Ｐｔ７１、Ｐｔ７２、Ｐｔ７３の経路長は、それぞれ２、１、２となる。
　図５Ｈに例示する経路パターンＰｎ０８は、４つの経路Ｐｔ８１－Ｐｔ８４を有する。経路Ｐｔ８１、Ｐｔ８２、Ｐｔ８３、Ｐｔ８４は、それぞれ位相決定済画素Ｏｐ０１、Ｏｐ１０、Ｏｐ０１、Ｏｐ２０を起点とする。経路Ｐｔ８１、Ｐｔ８２、Ｐｔ８３、Ｐｔ８４の経路長は、それぞれ１、１、３、２となる。
　図５Ｉに例示する経路パターンＰｎ０９は、４つの経路Ｐｔ９１－Ｐｔ９４を有する。経路Ｐｔ９１、Ｐｔ９２、Ｐｔ９３、Ｐｔ９４は、それぞれ位相決定済画素Ｏｐ０１、Ｏｐ０２、Ｏｐ１０、Ｏｐ１０を起点とする。経路Ｐｔ９１、Ｐｔ９２、Ｐｔ９３、Ｐｔ９４の経路長は、それぞれ１、３、１、３となる。
　図４、図５Ａ－図５Ｉは、終点のｘ座標とｙ座標の双方が、それぞれ起点のｘ座標とｙ座標以上となる経路のみを有し、起点の数が１以上３以下である経路パターンＰｎ００－Ｐｎ０９を例にしたが、これには限られない。終点のｘ座標が起点のｘ座標未満または終点のｙ座標が起点のｙ座標未満となる経路を有する経路パターンが設定されてもよいし、より起点の数が多い経路パターンが設定されてもよい。
　図６Ａ－図６Ｄ、図７Ａ－図７Ｃは、本実施形態に係る経路パターンの起点と終点のセットの例を示す。図６Ａ－図６Ｄ、図７Ａ－図７Ｃは、終点を中心とする３サンプル×３サンプルのｘ－ｙ平面内の領域に設定可能とする全てのセットを示す。但し、あるセットと起点の分布が終点を通るｘ軸またはｙ軸に対して対称なセットと、終点を中心として回転対称なセットについては、共通のセットとみなし、図示が省略されている。ここで、起点、終点は、それぞれ斜線、網掛けで示されている。また、経路の図示が省略されている。
　図６Ａ－図６Ｄは、起点の数がそれぞれ１－４となるセットを示す。起点の数が１、２、３、４個である場合、それぞれ２通り、４通り、８通り、１２通りのセットが設定可能である。
　図７Ａ－図７Ｃは、起点の数がそれぞれ５－７となるセットを示す。起点の数が５、６、７個である場合、それぞれ８通り、４通り、２通りのセットが設定可能である。
　なお、図４－図７Ｃは、起点は、終点からｘ方向とｙ方向の一方または双方（斜め方向）に隣接する位相決定済画像である場合である場合を例にしたが、終点から所定の範囲内であれば、これには限られない。例えば、起点は、終点から２サンプル以内の位相決定済画素であってもよい。
　また、図４、図５Ａ－図５Ｉに例示されるように、複数の積分経路がいずれも終点から所定の経路設定範囲（例えば、図８Ｂ、図８Ｃ、図９：経路設定範囲Ｐａ）内に含まれ、その一部または全部が所定の範囲外とならないように設定されることで、積分経路が過大となることが抑えられる。
　図４、図５Ａ－図５Ｉは、積分経路が、あるサンプル点からｘ方向またはｙ方向に隣接するサンプル点からなる素片（リンク）を、起点から終点に至るまで順次連接して形成される場合を例にしたが、これには限られない。個々の素片の起点は、その素片の終点から光学系の空間分解能の範囲内に配置されたサンプル点であればよい。素片は、例えば、ｘ方向とｙ方向の双方に交差する斜め方向に隣接する２個のサンプル点を結ぶ区間であってもよいし、互いに隣接せずに１以上のサンプル点を跨いだ（スキップ）２個のサンプル点を結ぶ区間であってもよい。また、積分経路の途中のサンプル点、つまり、起点と終点を除く積分経路上のサンプル点は、バルク位相エラーを既に定めたサンプル点であってもよいし、バルク位相エラーが未決定のサンプル点であってもよい。但し、位相勾配算出部１４０は、素片としてとりうる終点の位置ごとに、二次元領域内の各サンプル点について、そのサンプル点を素片の起点として終点との間のＯＣＴ信号の位相差を位相勾配として算出しておく。
（繰り返し演算）
　バルク位相エラー算出部１５０には、予め複数の経路パターン（例えば、経路パターンＰｎ００－Ｐｎ０９）を設定しておいてもよい。バルク位相エラー算出部１５０は、二次元領域内のサンプル点ごとに、そのサンプル点を終点として、起点として既にバルク位相エラーを定めたサンプル点との位置関係と対応関係を満たす経路パターンを、予め設定された経路パターンから選択する。バルク位相エラー算出部１５０は、選択した経路バターンで示される複数の経路のそれぞれを積分経路として上記の手法で終点におけるバルク位相エラーを算出する。そして、バルク位相エラー算出部１５０は、バルク位相エラーの算出する対象サンプル点を、バルク位相エラーを未決定の他のサンプル点に変更する。バルク位相エラー算出部１５０は、バルク位相エラーの算出と対象サンプル点の変更を、二次元領域内の全てのサンプル点についてバルク位相エラーが決定されるまで繰り返せばよい。対象サンプル点は、上記の位相決定対象画素に相当する。
　そこで、バルク位相エラー算出部１５０には、二次元領域内の１つのサンプル点を基準点とし、基準点に近接するサンプル点ほど優先されるようにバルク位相エラーを算出する演算順序（イテレーション）を予め定めておいてもよい。また、バルク位相エラー算出部１５０は、演算順序が第ｍ回（ｍは、１以上の整数）のサンプル点を終点とする対象サンプル点とするとき、演算順序が第ｍ－１回のサンプル点を起点として、その終点と起点が示す位置関係に対応する経路パターンを選択してもよい。ここで、基準点を、演算順序が第０回のサンプル点とし、そのバルク位相エラーを所定の基準値と設定しておく。また、同一の演算順序となるサンプル点は、必ずしも１個ではなく、２個以上となりうる。同一の演算順序となるサンプル点は、２個以上となりうるが、第ｍ回の２個以上のサンプル点の間では、バルク位相エラーの演算順序は、予め定めた順序でもよいし、ランダムに定めてもよい。
　図８Ａ－図８Ｃ、図９Ａに示す例では、バルク位相エラー算出部１５０には、二次元領域の原点（０，０）を基準点Ｏｐ００とし、原点からｘ方向とｙ方向に挟まれる斜め方向の距離の順に演算順序を設定しておく。同一の演算順序のサンプル点は、斜め方向に直交する線上に分布する。
　図８Ｂに示すように、演算順序が第１回のサンプル点は、（１，０）、（０，１）にそれぞれ所在するサンプル点が対象サンプル点Ｄｐ１０、Ｄｐ０１として設定される。このとき、バルク位相エラー算出部１５０は、対象サンプル点Ｄｐ１０、Ｄｐ０１のそれぞれについて、経路パターンＰｎ０１、Ｐｎ０２を選択することができる。但し、バルク位相エラー算出部１５０は、選択した経路パターンＰｎ０１、Ｐｎ０２のうち、経路Ｐｔ１１、Ｐｔ２１は、それぞれサンプル点が配列された二次元領域の範囲外を通過するため、無効な経路として棄却し、残りの有効な経路をバルク位相エラーの算出に用いる。経路別位相エラーを有効な経路間で合成する際に用いる重み係数は、有効な経路の数もしくは有効な経路の経路長に基づいて予め定めておいてもよい。
　図８Ｃに示すように、演算順序が第２回のサンプル点は、（１，１）、（２，０）、（１，０）にそれぞれ所在するサンプル点を対象サンプル点Ｄｐ１１、Ｄｐ２０、Ｄｐ０２として設定される。このとき、バルク位相エラー算出部１５０は、対象サンプル点Ｄｐ１０、Ｄｐ０１のそれぞれについて、経路パターンＰｎ０５、Ｐｎ０１、Ｐｎ０２を選択することができる。但し、バルク位相エラー算出部１５０は、選択した経路パターンＰｎ０１、Ｐｎ０２のうち、経路Ｐｔ１１、Ｐｔ２１は、それぞれサンプル点が配列された二次元領域の範囲外を通過するため、無効な経路として棄却し、残りの有効な経路をバルク位相エラーの算出に用いる。
　図９Ａに示すように、演算順序が第３回のサンプル点として、（２，１）、（１，２）にそれぞれ所在するサンプル点を対象サンプル点Ｄｐ２１、Ｄｐ１２として設定される。但し、図９Ａに示す例では、演算順序が第３回のサンプル点に、二次元領域において基準点Ｏｐ００から最も遠い対象サンプル点Ｄｐ２２も含まれる。このとき、バルク位相エラー算出部１５０は、対象サンプル点Ｄｐ２１、Ｄｐ１２、Ｄｐ２２のそれぞれについて、経路パターンＰｎ０５、Ｐｎ０５、Ｐｎ００を選択することができる。但し、バルク位相エラー算出部１５０は、対象サンプル点Ｄｐ２１、Ｄｐ１２、Ｄｐ２２のそれぞれについて選択した経路パターンＰｎ０５、Ｐｎ０５、Ｐｎ００のうち、経路Ｐｔ５４、Ｐｔ５３、Ｐｔ０３、Ｐｔ０４は、それぞれサンプル点が配列された二次元領域の範囲外を通過するため、無効な経路として棄却し、残りの有効な経路をバルク位相エラーの算出に用いる。
　二次元領域内のサンプル点間の演算順序は、図８Ａ－図８Ｃ、図９Ａに示す例には限られない。例えば、基準点とするサンプル点を二次元領域の中心とし、演算順序を基準点の周囲のサンプルのうち、基準点からの距離が小さいサンプルほど高くなるよう定めてもよい。そして、バルク位相エラー算出部１５０は、演算順序が第ｍ回（ｍは、１以上の整数）のサンプル点を終点とする対象サンプル点とするとき、演算順序が第ｍ－１回のサンプル点を起点として、その終点と起点が示す位置関係に対応する経路パターンを選択してもよい。
　図９Ｂは、一例として演算順序を示す経路Ｓｑが基準点を中心とする右回りの渦巻状に設定されているサンプル点Ａにおけるバルク位相エラーを算出する際、。バルク位相エラー算出部１５０は、サンプル点Ａを中心とする３×３点の領域を経路設定範囲ＰＡとして設定する。経路設定範囲ＰＡ内の位相決定済画素は基準点だけであるため、バルク位相エラー算出部１５０は、図５Ｂに例示される経路パターンを選択する。サンプル点Ｂにおけるバルク位相エラーを算出する際、バルク位相エラー算出部１５０は、サンプル点Ｂを中心とする３×３点の領域を経路設定範囲ＰＢとして設定する。経路設定範囲ＰＢ内の位相決定済画素は、サンプル点Ｂの真上、左上、左、左下の計４個であるため、バルク位相エラー算出部１５０は、図６Ｄの第１行左端に例示される経路パターンを選択する。
　このように、基準点からの距離が近接するほど優先されるようにサンプルごとの演算順序を定め、直前の演算順序を起点とする経路積分が実行されるので、二次元領域内において最後にバルク位相エラーを算出するまでの演算回数が減少する。そのため、二次元領域全体として、誤差の累積が抑制される。
　なお、対象サンプル点を終点とし、バルク位相エラーを決定済みのサンプル点を起点とする経路パターンは、その位置関係に応じて複数通り存在しうる。その場合には、バルク位相エラー算出部１５０には、複数通りの経路パターンのうちいずれか１通りを予め定めておいてもよいし、複数の通りの経路パターンから１通りをランダムに選択してもよい。
　二次元領域内のサンプル点間の演算順序は、これには限られない。例えば、バルク位相エラー算出部１５０は、各行におけるバルク位相エラーの算出処理（以下、ライン処理）として、位相決定済画素が二次元領域の終端のサンプル点に達するまで、位相決定済画素ごとに図５Ａに示す経路パターンを用いてバルク位相エラーを算出し、位相決定済画素をｘ軸方向に隣接するサンプル点に変更する処理を繰り返す。その後、バルク位相エラー算出部１５０は、図５Ｂに示す経路パターンを用いて、ｙ軸方向に隣接する次の行の始端のサンプル点におけるバルク位相エラーを算出する。そして、バルク位相エラー算出部１５０は、バルク位相の決定対象とする行が二次元領域の終端に達するまでライン処理を繰り返してもよい。
（位相の補正例）
　次に、算出されたバルク位相エラーによるＯＣＴ信号の位相の補正例について説明する。
　図１０Ａ－図１０Ｆは、本実施形態に係るＯＣＴ信号の位相の補正例を示す図である。
　図１０Ａは、ＯＣＴ信号に基づくＯＣＴ画像の一例を示す。図１０Ａに示すＯＣＴ画像は、観測可能領域の一部をなす観察対象平面の状態を示す。このＯＣＴ画像は、バルク位相エラーがまだ補償されていない生のＯＣＴ信号に基づく。
　図１０Ｄは、図１０Ａに係るＯＣＴ信号を空間周波数領域に変換して得られるスペクトルを示す。縦軸、横軸は、それぞれｘ方向、ｙ方向の空間周波数を示し、中心が原点を示す。スペクトルの強度は、濃淡で表されている。明るい部分ほど強度が高く、暗い部分ほど強度が低いことを示す。強度は原点に対して非対称であり、最も強度が高い周波数帯域の中心は原点からｘ方向に変位している。また、強度が高い部分の広がりが比較的大きい。
　図１０Ｂは、従来の位相差法を用いて算出されたバルク位相エラーの例を示す。位相は、濃淡で表されている。明るい部分ほど位相がπに近似し、暗い部分ほど位相が－πに近似することを示す。バルク位相エラーは、ｘ方向に対しては比較的規則的に変化するが、ｙ方向の変化がｘ方向よりも不規則となる。
　図１０Ｅは、図１０Ｂに示すバルク位相エラーを用いて補正されたスペクトルを示す。スペクトルの強度はｘ方向に対しては、原点を通るｙ軸に対して対称となり、ｘ方向の広がりが緩和する。しかしながら、スペクトル強度のｙ方向の広がりが却って拡大している。このことは、ｙ方向に対する補正が不適切なことを示す。
　図１０Ｃは、本実施形態に係る信号処理装置１００が算出したバルク位相エラーの例を示す。位相は、濃淡で表されている。このバルク位相エラーは、図８Ａ－図８Ｃ、図９を用いて説明した手法を用いて算出されたものである。バルク位相エラーは、ｘ方向、ｙ方向の両者に対して規則的に変化し、比較的滑らかな縞模様を示す。
　図１０Ｆは、図１０Ｃに示すバルク位相エラーを用いて補正されたＯＣＴ信号のスペクトルを示す。
　信号補正部１６０は、式（６）を用いてＯＣＴ信号を補正することにより、ＯＣＴ信号に含まれるバルク位相エラーを除去することができる。式（６）は、三次元サンプル点（ｘ，ｙ，ｚ）における信号値のスペクトルＡ（ｘ，ｙ，ｚ）の位相からバルク位相エラーφｂ （ｔ）（ｘ，ｙ）を差し引いて、補正後の信号値のスペクトルＡ’（ｘ，ｙ，ｚ）を算出することを示す。スペクトルの強度はｘ方向、ｙ方向ともに原点に対して対称となり、広がりが緩和する。このことは、本実施形態により従来の位相差法よりも位相エラーに生じる方向依存性を緩和し、より適切なスペクトルの補正を実現できることを示す。
Figure JPOXMLDOC01-appb-M000006
　次に、試料ごとのＯＣＴ信号の別の補正例について説明する。図１１Ａ、図１１Ｂ、図１１Ｃ、図１１Ｄは、４種類の試料に対して観測されたｘ方向の空間周波数ｆｘの平均値、ｙ方向の空間周波数ｆｙの平均値、空間周波数ｆｘの標準偏差、空間周波数ｆｙの標準偏差の例をそれぞれ示す。これらの平均値、標準偏差は、補正前の原データ、単純積分（simple integration）により算出されたバルク位相エラーで補正した場合、本実施形態における提案法により算出されたバルク位相エラーで補正した場合のそれぞれについて算出した。Ｓｍ０１、Ｓｍ０３は、それぞれニワトリの胸筋（chicken breast muscle）の異なる部位を試料とした場合を示す。Ｓｍ０２、Ｓｍ０４は、それぞれブタの心臓組織（porcine heart tissue）の異なる部位を試料とした場合を示す。Ａｖｇは、Ｓｍ０１－Ｓｍ０４間の平均値を示す。
　図１１Ａ、図１１Ｂは、原データについては、試料間で平均値が著しくばらつき、単純積分、提案法について、平均値はほぼ０となり、原データよりも有意に低減することを示す。空間周波数ｆｙの平均値については、提案法よりも単純積分の方がむしろ、ばらつきが抑えられる。
　他方、図１１Ｃは単純積分により空間周波数ｆｘの標準偏差が原データよりも少なくなるが、図１１Ｄは単純積分により空間周波数ｆｙの標準偏差が原データよりも増加する。このことは、ｘ方向よりもｙ方向にスペクトルが広がることを示す。これに対し、提案法によれば原データとの間で空間周波数ｆｘ、ｆｙの標準偏差の有意な変化が認められない。このことは、提案法によれば特定方向へのスペクトルの偏りをもたらさずにバルク位相エラーを補正できることを示す。
　なお、信号補正部１６０は、補正後のＯＣＴ信号に対して、再焦点化してもよい（ディジタルリフォーカス、digital refocusing）。再焦点化において、信号補正部１６０は、次の各ステップの処理を実行する。
（ステップＳ１１）信号補正部１６０は、補正後の二次元領域のＯＣＴ信号（例えば、正面ＯＣＴ信号、en face OCT slice）に対して二次元離散フーリエ変換（ＤＦＴ：Discrete Fourier Transform）を行って、空間周波数領域における複素スペクトル（complex spatial frequency spectrum）を算出する。
（ステップＳ１２）信号補正部１６０は、複素スペクトルを構成する空間周波数（ｆｘ，ｆｙ）ごとの変換係数に、その空間周波数に対する位相フィルタＨ－１（ｆｘ，ｆｙ）を乗算することにより、乗算により得られる乗算値を含むフィルタ処理後のスペクトルを算出する。位相フィルタＨ－１（ｆｘ，ｆｙ）の例を、式（７）に示す。
Figure JPOXMLDOC01-appb-M000007
　式（７）において、ｆｒは、動径方向の空間周波数を示す。ｆｒは、ｘ方向の空間ｆｘとｙ方向の空間周波数ｆｙの二乗和の平方根に相当する。λｃは、入射光の中心波長を示す。ｚ０は、焦点ずれ量（amount of defocus）を示す。焦点ずれ量ｚ０は、一般に光学系に依存する。位相フィルタＨ－１（ｆｘ，ｆｙ）は、空間周波数が高いほど顕著に表れる位相のずれを補償することで、焦点ずれを緩和または除去することを可能にする。かかる位相のずれは、典型的には、「ぼけ」としてＯＣＴ画像に表れる。
（ステップＳ１３）信号補正部１６０は、フィルタ処理後のスペクトルに対して二次元逆フーリエ変換（ＩＤＦＴ：Inverse Discrete Fourier Transform）を行ってフィルタ処理後のＯＣＴ信号を再焦点化したＯＣＴ信号として生成する。そして、信号補正部１６０は、再焦点化したＯＣＴ信号を画像処理部１７０に出力する。画像処理部１７０によれば、信号補正部１６０から入力されるＯＣＴ信号に基づき、より鮮明なＯＣＴ画像を取得することができる。
　以上に説明したように、本実施形態に係る信号処理装置１００は、位相勾配算出部１４０とバルク位相エラー算出部１５０を備える。位相勾配算出部１４０は、試料の状態を表すＯＣＴ信号の位相勾配を所定の平面（例えば、二次元領域）に配列されたサンプル点ごとに算出する。バルク位相エラー算出部１５０は、バルク位相エラーが決定済のサンプル点である起点からバルク位相エラーが未決定のサンプル点である終点までの複数の経路のそれぞれについて、サンプル点ごとの位相勾配をそれぞれの経路に沿って積分して前記終点における経路別位相エラーを算出し、経路別位相エラーを複数の経路間で合成して終点におけるバルク位相エラーを定める。バルク位相エラー算出部１５０は、バルク位相エラーとして、例えば、経路別位相エラーの経路間の平均値を前記バルク位相エラーとして算出する。
　この構成によれば、それぞれ経路によって異なりうる誤差のランダムな成分が、経路間で合成されることで相対的に低減する。そのため、経路によらない成分がバルク位相エラーとして、単純に位相勾配を累積するよりも正確に推定される。
　また、本実施形態において、バルク位相エラー算出部１５０は、経路間の平均値として、複数の経路それぞれの経路長に基づく重み係数を用いた経路別位相エラーの加重平均値を算出してもよい。
　この構成によれば、経路長に依存した誤差の影響を考慮して、より正確にバルク位相エラーを算出することができる。例えば、経路長が大きい経路ほど重み係数を小さくすることで、経路長が大きいほど誤差の影響が強い場合、かかる経路の経路別位相エラーの寄与を小さくすることができる。逆に、経路長が大きい経路ほど重み係数を大きくすることで、経路長が大きいほど誤差の影響が弱い場合、かかる経路の経路別位相エラーの寄与を大きくすることができる。
　また、本実施形態において、バルク位相エラーの算出に用いられる複数の経路のそれぞれは、終点から所定の範囲内に含まれてもよい。
　この構成によれば、複数の経路のそれぞれが所定の範囲内に限定されるため、バルク位相エラーの算出に係る経路数が過大になることが回避される。
　また、本実施形態において、バルク位相エラーの算出に用いられる複数の経路のそれぞれは、第ｎ（ｎは、１以上かつＮ－１以下の整数）サンプル点から第ｎ＋１サンプル点の区間である第ｎ素片（Ｎは、２以上の整数）を含み、第ｎ＋１サンプル点は、第ｎサンプル点からＯＣＴ信号を取得した光学系の空間分解能以下の範囲内のサンプル点であってもよい。
　これにより、複数の経路のそれぞれは、ＯＣＴ信号の信号値に有意差が生じない２つのサンプル点からなる１個または複数の素片を連接して形成されるため、算出されるバルク位相エラーからＯＣＴ信号の信号成分による寄与を抑制することができる。
　また、本実施形態において、バルク位相エラー算出部１５０には、起点と終点の位置関係に対応付けて複数の経路を示す経路パターンを予め複数通り設定しておき、バルク位相エラーの算出対象とする対象サンプル点を終点として、バルク位相エラーを決定済みのサンプル点を起点とする位置関係を満たす経路パターンで示される複数の経路を選択してもよい。
　これにより、二次元領域内で対象サンプル点ごとに、バルク位相エラーを算出する際に用いる複数の経路を、バルク位相エラーの算出済のサンプル点との関係で、その都度探索するよりも効率的に選択することができる。
　また、本実施形態において、バルク位相エラー算出部１５０には、二次元平面（例えば、二次元領域）における１つのサンプル点を基準点とし、基準点に近接するサンプル点ほど優先されるように、バルク位相エラーを算出する順序を予め定めておき、定めた順序が第ｍ回（ｍは、１以上の整数）のサンプル点を前記対象サンプル点とするとき、定めた順序が第ｍ－１回のサンプル点を起点とする位置関係を満たす経路パターンで示される複数の経路を選択してもよい。
　これにより、順次算出されるバルク位相エラーに累積される誤差の影響が、演算順序を任意とする場合よりも抑制されるため、二次元平面全体として、バルク位相エラーをより高い精度で算出することができる。
　また、本実施形態に係る信号処理装置１００は、信号補正部１６０を備えてもよい。信号補正部１６０は、終点におけるバルク位相エラーを、終点と前記平面に射影された位置が共通する三次元空間内に配列されたサンプル点におけるＯＣＴ信号の信号値の位相から減算してもよい。
　これにより、三次元のＯＣＴ信号に含まれるバルク位相エラーが補償されるため、バルク位相エラーによるスペクトルの偏りや広がりを低減し、より高品質の光干渉断層画像を得ることができる。
　以上、図面を参照してこの発明の実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。
　例えば、上記の説明では、信号処理装置１００が光干渉断層計１の一部である場合を例にしたが、これには限られない。信号処理装置１００は、光干渉断層計１から独立し、光学系を備えない単一の機器であってもよい。その場合、信号処理装置１００の制御部１１０において光学系制御部１２０が省略されてもよい。検出信号取得部１３０は、光学系に限られず、データ蓄積装置、ＰＣなど、他の機器から有線または無線で、例えば、ネットワークを経由して検出信号またはＯＣＴ信号を取得してもよい。
　また、信号処理装置１００は、上記の操作入力部と表示部を備えてもよいし、それらの一方または両方が省略されてもよい。
　信号処理装置１００の制御部１１０において画像処理部１７０が省略されてもよい。その場合、信号処理装置１００の制御部１１０は、信号補正部１６０が生成した補正ＯＣＴ信号を、データ蓄積装置、ＰＣ、画像処理装置など、他の機器に出力してもよい。出力先とする機器は、画像処理部１７０と同様の機能、つまり、信号処理装置１００から入力される補正ＯＣＴ信号に基づいて画像データを生成し、生成した画像データに基づくＯＣＴ画像を表示する機能を有してもよい。
　上記の例では、位相勾配もしくはバルク位相エラーの算出に係るサンプル点が配列された二次元領域がｘ－ｙ平面に設定される場合を例にしたが、これには限られない。観測可能領域をなす三次元空間に配列された三次元サンプル点と、その三次元サンプル点を射影して得られる二次元サンプル点との対応関係を形成できる二次元平面に設定されればよい。かかる二次元平面は、例えば、ｙ－ｚ平面、ｚ－ｘ平面などであってもよい。かかる二次元平面がｙ－ｚ平面である場合には、三次元サンプル点（ｘ，ｙ，ｚ）は、二次元サンプル点（ｙ，ｚ）に対応付けられる。また、二次元領域として、ＯＣＴ画像を構成する画素に対応するサンプル点が配列された観察対象平面が用いられてもよい。信号補正部１６０は、バルク位相エラーに基づくＯＣＴ信号の位相の補正を、個々の観察対象平面に対して実行してもよい。
　なお、上述した実施形態における信号処理装置１００の一部、例えば、制御部１１０の全部または一部をコンピュータで実現するようにしてもよい。その場合、この制御機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現してもよい。なお、ここでいう「コンピュータシステム」とは、信号処理装置１００に内蔵されたコンピュータシステムであって、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ－ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでもよい。また上記プログラムは、前述した機能の一部を実現するためのものであってもよく、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであってもよい。
　また、上述した実施形態及び変形例における信号処理装置１００の一部、または全部を、ＬＳＩ（Large Scale Integration）等の集積回路として実現してもよい。信号処理装置１００の各機能ブロックは個別にプロセッサ化してもよいし、一部、または全部を集積してプロセッサ化してもよい。また、集積回路化の手法はＬＳＩに限らず専用回路、または汎用プロセッサで実現してもよい。また、半導体技術の進歩によりＬＳＩに代替する集積回路化の技術が出現した場合、当該技術による集積回路を用いてもよい。
１…光干渉断層計、１０…光源、２０…ビームスプリッタ、３０ａ、３０ｂ、５０ａ、５０ｂ…コリメータ、４０…参照鏡、６０ａ、６０ｂ…ガルバノミラー、７０…分光器、１００…信号処理装置、１２０…光学系制御部、１３０…検出信号取得部、１４０…位相勾配算出部、１５０…バルク位相エラー算出部、１６０…信号補正部、１７０…画像処理部

Claims (10)

  1. 　試料の状態を表す光干渉断層信号の光の照射方向に交差する平面における位相勾配を前記平面に配列されたサンプル点ごとに算出する位相勾配算出部と、
    　バルク位相エラーが決定済のサンプル点である起点からバルク位相エラーが未決定のサンプル点である終点までの複数の経路のそれぞれについて、サンプル点ごとの前記位相勾配をそれぞれの経路に沿って積分して前記終点における経路別位相エラーを算出し、
    　前記経路別位相エラーを前記複数の経路間で合成して前記終点におけるバルク位相エラーを定めるバルク位相エラー算出部と、を備える
    　信号処理装置。
  2. 　前記バルク位相エラー算出部は、
    　前記経路別位相エラーの経路間の平均値を前記バルク位相エラーとして算出する
    　請求項１に記載の信号処理装置。
  3. 　前記バルク位相エラー算出部は、
    　前記平均値として、前記複数の経路それぞれの経路長に基づく重み係数を用いた前記経路別位相エラーの加重平均値を算出する
    　請求項２に記載の信号処理装置。
  4. 　前記複数の経路のそれぞれは、前記終点から所定の範囲内に含まれる
    　請求項１から請求項３のいずれか一項に記載の信号処理装置。
  5. 　前記複数の経路のそれぞれは、第ｎ（ｎは、１以上かつＮ－１以下の整数）サンプル点から第ｎ＋１サンプル点の区間である第ｎ素片（Ｎは、２以上の整数）を含み、
    　前記第ｎ＋１サンプル点は、前記第ｎサンプル点から前記光干渉断層信号を取得した光学系の空間分解能以下の範囲内のサンプル点である
    　請求項１から請求項４のいずれか一項に記載の信号処理装置。
  6. 　前記バルク位相エラー算出部には、
    　前記起点と前記終点の位置関係に対応付けて前記複数の経路を示す経路パターンを予め複数通り設定しておき、
    　バルク位相エラーの算出対象とする対象サンプル点を終点として、前記位置関係を満たす経路パターンで示される複数の経路を選択する
    　請求項１から請求項５のいずれか一項に記載の信号処理装置。
  7. 　前記バルク位相エラー算出部には、
    　前記平面における１つのサンプル点を基準点とし、
    　前記基準点に近接するサンプル点ほど優先されるように、バルク位相エラーを算出する順序を予め定めておき、
    　前記順序が第ｍ回（ｍは、１以上の整数）のサンプル点を前記対象サンプル点とするとき、前記順序が第ｍ－１回のサンプル点を起点として、前記位置関係を満たす経路パターンで示される前記複数の経路を選択する
    　請求項６に記載の信号処理装置。
  8. 　前記終点におけるバルク位相エラーを、前記終点と前記平面に射影された位置が共通する三次元空間内に配列されたサンプル点における前記光干渉断層信号の信号値の位相から減算する信号補正部を備える
    　請求項１から請求項７のいずれか一項に記載の信号処理装置。
  9. 　信号処理装置における方法であって、
    　試料の状態を表す光干渉断層信号の光の照射方向に交差する平面における位相勾配を前記平面に配列されたサンプル点ごとに算出する位相勾配算出過程と、
    　バルク位相エラーが決定済のサンプル点である起点からバルク位相エラーが未決定のサンプル点である終点までの複数の経路のそれぞれについて、サンプル点ごとの前記位相勾配をそれぞれの経路に沿って積分して前記終点における経路別位相エラーを算出し、
    　前記経路別位相エラーを前記複数の経路間で合成して前記終点におけるバルク位相エラーを定めるバルク位相エラー算出過程と、を有する
    　信号処理方法。
  10. 　信号処理装置のコンピュータに、
    　試料の状態を表す光干渉断層信号の光の照射方向に交差する平面における位相勾配を前記平面に配列されたサンプル点ごとに算出する位相勾配算出手順と、
    　バルク位相エラーが決定済のサンプル点である起点からバルク位相エラーが未決定のサンプル点である終点までの複数の経路のそれぞれについて、サンプル点ごとの前記位相勾配をそれぞれの経路に沿って積分して前記終点における経路別位相エラーを算出し、
    　前記経路別位相エラーを前記複数の経路間で合成して前記終点におけるバルク位相エラーを定めるバルク位相エラー算出手順と、
    　を実行させるための信号処理プログラム。
PCT/JP2020/042459 2019-11-14 2020-11-13 信号処理装置、信号処理方法および信号処理プログラム WO2021095852A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
US17/776,178 US20220390368A1 (en) 2019-11-14 2020-11-13 Signal processing device, signal processing method, and signal processing program
JP2021556177A JPWO2021095852A1 (ja) 2019-11-14 2020-11-13

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2019-206506 2019-11-14
JP2019206506 2019-11-14

Publications (1)

Publication Number Publication Date
WO2021095852A1 true WO2021095852A1 (ja) 2021-05-20

Family

ID=75913004

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2020/042459 WO2021095852A1 (ja) 2019-11-14 2020-11-13 信号処理装置、信号処理方法および信号処理プログラム

Country Status (3)

Country Link
US (1) US20220390368A1 (ja)
JP (1) JPWO2021095852A1 (ja)
WO (1) WO2021095852A1 (ja)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002504671A (ja) * 1998-02-23 2002-02-12 ザイゴ コーポレイション 空気の屈折率および光路長の影響を測定するための干渉測定器および方法
JP2013140077A (ja) * 2012-01-05 2013-07-18 Univ Of Tsukuba 波長走査型光干渉断層計及びその位相安定化プログラム
US20140375792A1 (en) * 2011-12-09 2014-12-25 Massachusetts Institute Of Technology Systems and methods for self-referenced quantitative phase microscopy
JP2016504947A (ja) * 2013-02-01 2016-02-18 カール ツアイス メディテック アクチエンゲゼルシャフト 干渉撮像におけるサブアパーチャベースの収差測定及び補正用のシステム及び方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002504671A (ja) * 1998-02-23 2002-02-12 ザイゴ コーポレイション 空気の屈折率および光路長の影響を測定するための干渉測定器および方法
US20140375792A1 (en) * 2011-12-09 2014-12-25 Massachusetts Institute Of Technology Systems and methods for self-referenced quantitative phase microscopy
JP2013140077A (ja) * 2012-01-05 2013-07-18 Univ Of Tsukuba 波長走査型光干渉断層計及びその位相安定化プログラム
JP2016504947A (ja) * 2013-02-01 2016-02-18 カール ツアイス メディテック アクチエンゲゼルシャフト 干渉撮像におけるサブアパーチャベースの収差測定及び補正用のシステム及び方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
FINGLER, J ET AL.: "Volumetric microvascular imaging of human retina using optical coherence tomography with a novel motion contrast technique", OPTICS EXPRESS, vol. 17, no. 24, 23 November 2009 (2009-11-23), XP055291503, DOI: 10.1364/OE.17.022190 *
OIKAWA, K ET AL.: "Bulk-phase-error correction for phase-sensitive signal processing of optical coherence tomography", BIOMEDICAL OPTICS EXPRESS, vol. 11, no. 10, 25 September 2020 (2020-09-25), pages 5886, XP055823385, DOI: https://doi.org/10.1364/B0E.396666 *
SHEMONSKI, N. D ET AL.: "Three-dimensional motion correction using speckle and phase for in vivo computed optical interferometric tomography", BIOMEDICAL OPTICS EXPRESS, vol. 5, no. 12, 1 December 2014 (2014-12-01), XP055823380, DOI: https:// doi.org/10.1364/BOE.5.004131 *
SPAHR, H ET AL.: "Interferometric detection of 3D motion using computational subapertures in optical coherence tomography", OPTICS EXPRESS, vol. 26, no. 15, 9 July 2018 (2018-07-09), pages 18803, XP055823381, DOI: https://doi.org/10.1364/0E.26.018803 *
YASUNO, Y ET AL.: "Non-iterative numerical method for laterally superresolving Fourier domain optical coherence tomography", OPTICS EXPRESS, vol. 14, no. 3, 6 February 2006 (2006-02-06), pages 1006, XP055823384 *

Also Published As

Publication number Publication date
JPWO2021095852A1 (ja) 2021-05-20
US20220390368A1 (en) 2022-12-08

Similar Documents

Publication Publication Date Title
JP4389032B2 (ja) 光コヒーレンストモグラフィーの画像処理装置
US9354038B2 (en) Swept source optical coherence tomography and method for stabilizing phase thereof
US7301644B2 (en) Enhanced optical coherence tomography for anatomical mapping
EP2149776B1 (en) Optical coherence tomographic imaging method and optical coherence tomographic imaging apparatus
JP4461259B2 (ja) 光断層画像の処理方法
JP6360065B2 (ja) スペクトル領域干渉法における信号処理方法および装置、並びにスペクトル領域光コヒーレンストモグラフィの方法および装置
US9615736B2 (en) Optical interference tomographic apparatus, and method for controlling optical interference tomographic apparatus
JP6798095B2 (ja) 光コヒーレンストモグラフィ装置、及びそれに用いる制御プログラム
JP6840520B2 (ja) 画像処理装置、撮像装置、画像処理方法及びプログラム
US9243887B2 (en) Lateral distortion corrected optical coherence tomography system
US20200129068A1 (en) Intraoral oct with color texture
US10524663B2 (en) Phase measurement, analysis, and correction methods for coherent imaging systems
US20180000341A1 (en) Tomographic imaging apparatus, tomographic imaging method, image processing apparatus, image processing method, and program
JP6431559B2 (ja) 光干渉断層撮影におけるモーションアーチファクトの除去のための方法およびシステム
US20220133446A1 (en) High-speed, dental optical coherence tomography system
US10760893B2 (en) Apparatus and method for processing the signal in master slave interferometry and apparatus and method for master slave optical coherence tomography with any number of sampled depths
Ksenofontov et al. Numerical method for axial motion artifact correction in retinal spectral-domain optical coherence tomography
WO2021095852A1 (ja) 信号処理装置、信号処理方法および信号処理プログラム
JP6214020B2 (ja) 光断層イメージング法、その装置およびプログラム
WO2016098850A1 (ja) 画像処理装置、レーザ照射システム、画像処理方法、および画像処理プログラム
Gelikonov et al. Numerical method for axial motion correction in optical coherence tomography
JP2018140004A (ja) 撮像装置、撮像方法およびプログラム
WO2021095826A1 (ja) 画像処理装置、画像処理方法およびプログラム
JP6760310B2 (ja) 画像データ処理装置、画像データ処理方法、および画像データ処理プログラム
JP5746741B2 (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: 20887825

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2021556177

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20887825

Country of ref document: EP

Kind code of ref document: A1