WO2016031697A1 - Cartilage diagnosis device and diagnostic probe - Google Patents

Abstract

　A cartilage diagnosis device 1 is provided with: an optical unit 2 including an OCT optical system; a probe including an optical mechanism connected to the optical unit 2 and configured so that a distal-end part thereof is insertable into an articular cavity, the optical mechanism for guiding light from the optical unit 2 to cartilage and causing the light to scan, and a weight mechanism capable of loading a predetermined weight on the cartilage; a control computation unit 6 for controlling driving of the weight mechanism and the optical mechanism, processing an optical interference signal outputted from the optical unit 2 on the basis of the driving of the weight mechanism and the optical mechanism, performing computation whereby a change in a predetermined dynamic characteristic quantity accompanying deformation in the cartilage due to weight loading is correlated with a tomographic position of the cartilage, and computing the degree of damage of cartilage tissue on the basis of the change in the dynamic characteristic quantity; and a display device 8 for displaying the degree of damage of the cartilage tissue in the form of a tomographic visualization.

Description

Cartilage diagnostic apparatus and diagnostic probe

The present invention relates to an apparatus for diagnosing the degree of cartilage degeneration.

Cartilage plays important roles such as relaxation of load impact and improvement of joint slidability, but it is a tissue that has no blood circulation and is difficult to self-heal. Many elderly people develop osteoarthritis (Osteoarthritis: hereinafter referred to as “OA”) due to cartilage wear, and the establishment of a diagnostic treatment method is required especially in countries with an aging society. . Articular cartilage contains 80% moisture and 20% matrix, which is composed of collagen and proteoglycans. In particular, it is believed that proteoglycans constrained by the collagen fibers determine the flow characteristics inside the cartilage and are largely involved in the excellent viscoelastic properties of the cartilage. OA is caused by a loss of the viscoelastic properties of the cartilage.

With the recent development of medical diagnostic technology, optical coherence tomography (Optical Coherence Tomography: hereinafter referred to as “OCT”) has been developed. According to this OCT, a microtomography can be visualized inside a living tissue in a non-invasive and non-contact manner. Further, the acquisition rate of the two-dimensional OCT tomographic image is equal to or higher than the video rate and has a high temporal resolution. In view of this, a technique has been proposed in which the OCT is used to visualize the tomographic characteristics of cartilage dynamics (Non-Patent Document 1). According to this technique, there is a possibility that normal cartilage and degenerated cartilage can be distinguished by applying a predetermined stress to the cartilage and visualizing the stress relaxation process.

Although the technology described in Non-Patent Document 1 provides a powerful foothold for micromechanics visualization diagnosis related to cartilage, there is room for improvement in terms of diagnostic efficiency and diagnostic accuracy. Moreover, it did not disclose a specific device configuration in the case where the cartilage of a living human body is a diagnosis target. For this reason, it has not yet reached practical use.

The present invention has been made in view of such problems, and one of its purposes is to make cartilage diagnosis using OCT more practical.

A certain aspect of the present invention is a cartilage diagnostic apparatus for diagnosing articular cartilage. This cartilage diagnostic apparatus is configured to include an optical unit including an optical system that uses optical coherence tomography, and is connected to the optical unit, and a tip portion thereof is configured to be inserted into a joint cavity, and guides light from the optical unit to the cartilage for scanning. And a probe including a load device for applying a predetermined deformation energy to the cartilage, and controls the drive of the load device and the optical mechanism, and is output from the optical unit based on the drive. The optical interference signal is processed, and a change in the predetermined mechanical feature amount due to the deformation inside the cartilage due to the application of deformation energy is calculated in association with the tomographic position of the cartilage. Based on the change in the mechanical feature amount, the cartilage tissue A control calculation unit that calculates the degree of damage of the cartilage, and a display device that displays the degree of damage of the cartilage tissue in a manner of visualizing a tomogram.

According to this aspect, a diagnosis target is obtained by a calculation process using the fact that there is a correspondence between the change in the mechanical feature amount generated in the cartilage when a predetermined deformation energy is applied and the degree of damage to the cartilage tissue. The degree of cartilage damage is visualized by tomography. For this reason, a doctor or the like can easily perform cartilage diagnosis by confirming the display of the tomographic visualization. That is, cartilage diagnosis using OCT can be more practically used.

Another aspect of the present invention is a diagnostic probe. This diagnostic probe is a probe that enables diagnosis of a target site by being connected to an optical unit including an optical system that uses optical coherence tomography, and includes a tip portion configured to be able to contact the target site, An optical mechanism for guiding the light from the optical unit to the cartilage to scan, and a load device for applying a predetermined deformation energy to the target site.

By bringing the diagnostic probe of this aspect into contact with the cartilage as the target site and operating the optical unit, it becomes possible to visualize the degree of damage of the cartilage described above. That is, cartilage diagnosis using OCT can be more practically used.

According to the present invention, cartilage diagnosis using OCT can be more practically used.

It is a figure which represents roughly the structure of the cartilage diagnostic apparatus which concerns on 1st Example. It is a figure which represents roughly the structure of the diagnostic probe which comprises a cartilage diagnostic apparatus. It is a figure which shows roughly the process sequence by a recursive cross correlation method. It is a figure which shows the process sequence by a subpixel analysis roughly. It is a figure which shows the result obtained by each process. It is a figure showing the process sequence of a cartilage diagnostic process. It is a figure showing the process sequence of a cartilage diagnostic process. It is a figure showing the composition of an experimental device roughly. It is a figure showing the test piece used for experiment. It is a graph which shows the load state applied to the test piece in the experiment. It is the figure which visualized the tomographic change of the strain rate distribution by the stress relaxation test. It is the figure which visualized the tomographic change of the strain rate distribution by the stress relaxation test. It is a figure showing the comparison of the attenuation coefficient of a normal cartilage and a simulation degenerated cartilage. It is a flowchart which shows the flow of the cartilage diagnostic process performed by a control calculating part. It is a flowchart which shows in detail the modification | denaturation degree calculation presentation process of S36 in FIG. It is a flowchart which shows the modification | denaturation degree calculation presentation process which concerns on 2nd Example in detail.

One embodiment of the present invention is a cartilage diagnostic apparatus for diagnosing articular cartilage. This cartilage diagnostic apparatus takes a tomographic image of cartilage using OCT while applying a predetermined deformation energy (load) to the cartilage, and calculates the degree of damage (degeneration degree) from the behavior of the cartilage tissue with respect to the load. Is. Since the calculation result is visualized on the display device in the form of a tomographic image of cartilage, a doctor or the like can make a cartilage diagnosis by viewing the tomographic image.

This cartilage diagnostic apparatus includes an optical unit including an OCT optical system and a probe connected to the optical unit. This probe has a tip portion that can be inserted into a joint cavity, and includes an optical mechanism and a load mechanism (functioning as a “load device”) for obtaining a physical quantity used for OCT calculation processing. The tip of the probe may be inserted into the joint cavity, for example, through a syringe needle. A predetermined load is applied to the cartilage by driving the load mechanism in a state where the tip of the probe is in contact with the surface (or the vicinity thereof) of the cartilage. The degree of cartilage degeneration can be evaluated by tomographically measuring the response of the cartilage tissue to the load by OCT. This is because, as will be described later, the response appears differently depending on the progress of the degree of denaturation. In order to acquire the behavior of the cartilage tissue as a tomographic image, the optical mechanism irradiates and scans the cartilage with object light from the optical unit, acquires the reflected light, and sends it to the optical unit.

The optical unit combines the reflected light and the reference light, and sends the optical interference signal to the control calculation unit. The control calculation unit processes the optical interference signal while controlling the driving of the load mechanism and the optical mechanism, and associates the change in the mechanical feature amount due to the deformation in the cartilage due to the load load described above with the tomographic position of the cartilage. To calculate. The “mechanical feature amount” here may be obtained based on the spatial distribution of the deformation vector of the cartilage tissue. For example, a deformation rate vector obtained by temporally differentiating the deformation vector or a strain rate tensor obtained by further spatially differentiating the deformation rate vector may be used.

The load loading method using a load mechanism may be based on a stress relaxation method in which a constant strain is applied to a measurement target (cartilage) and a time change of stress is measured. Or based on the dynamic viscoelasticity method which gives the dynamic strain with respect to a measuring object, and measures the maximum value and phase difference of stress. Alternatively, it may be based on a creep method in which a time-dependent change in strain is measured by applying a certain amount of stress to the measurement object. The control calculation unit calculates the degree of damage of the cartilage tissue based on the change in the mechanical feature value, and visualizes it on the screen of the display device.

The control calculation unit may execute the following processing based on the stress relaxation method. That is, the stress is relaxed after a predetermined load is applied to the cartilage by driving the load mechanism, and the deformation velocity vector corresponding to the tomographic position of the cartilage is determined based on the tomographic image data before and after the stress is taken by the optical unit. You may calculate as a mechanical feature-value. Then, the tomographic distribution of the attenuation coefficient obtained from the change of the deformation velocity vector may be visualized as a tomographic equivalent as the distribution of the damage degree of the cartilage tissue. The “front and back tomographic image data” here may be based on two tomographic images that are continuously captured, or may be based on two tomographic images that are not consecutive but are captured at a predetermined time.

Alternatively, the control calculation unit calculates a strain rate tensor by spatially differentiating the deformation rate vector calculated as a mechanical feature amount, and the fault distribution of the strain rate attenuation coefficient is equivalent to the distribution of the degree of damage of the cartilage tissue. It may be made visible. By calculating the strain rate tensor in this way, physical implications such as compression, expansion and shearing of the cartilage tissue can be given to the fault distribution of the attenuation coefficient. From this physical meaning, the degree of cartilage degeneration can be directly transmitted to a doctor or the like as a mechanical deformability.

The control calculation unit may calculate the deformation velocity vector with sub-pixel accuracy based on the cross-correlation of the tomographic image data acquired sequentially, the displacement of the luminance gradient, and the deformation of the luminance pattern. Since the resolution of OCT is limited, by introducing such sub-pixel analysis, cartilage diagnosis with higher accuracy can be realized.

The control calculation unit may execute the following processing based on the dynamic viscoelasticity method. That is, a predetermined dynamic load (dynamic force) is applied to the cartilage by driving the load mechanism, and each cartilage of the cartilage is determined based on the tomographic image data taken by the optical unit during the vibration of the cartilage due to the dynamic load. The deformation velocity vector at the tomographic position is calculated as a mechanical feature, and the viscoelastic tomographic distribution of cartilage is calculated based on at least one of the amplitude and phase difference of the specific deformation velocity vector that is the excitation frequency component of the deformation velocity vector. Alternatively, the viscoelastic fault distribution may be visualized as equivalent to the distribution of the damage degree of the cartilage tissue. This is because the viscoelasticity of the cartilage changes depending on the degree of degeneration of the cartilage. The “front and back tomographic image data” here may be based on two tomographic images that are continuously captured, or may be based on two tomographic images that are not consecutive but are captured at a predetermined time.

Alternatively, the strain rate tensor at each tomographic position of the cartilage is calculated as a mechanical feature, and the viscoelasticity of the cartilage is determined based on at least one of the amplitude and phase difference of the specific strain rate tensor that is the excitation frequency component of the strain rate tensor. The fault distribution may be calculated.

Further, the load mechanism may include a sensor capable of measuring a load (force) applied to the cartilage. The control calculation unit calculates the tomographic distribution of the viscoelasticity of the cartilage based on at least one of the amplitude and phase difference of the excitation frequency component of the calculated mechanical feature and the load value measured by the sensor, The tomographic distribution may be visualized as equivalent to the distribution of the degree of damage of the cartilage tissue.

The probe can also be applied to diagnosis of a target site other than cartilage, for example, for diagnosing internal organs during surgery. That is, the probe may be configured as a diagnostic probe that enables diagnosis of a predetermined target site by being connected to an optical unit including an OCT optical system. This probe has a tip that is configured to be able to contact the target part, an optical mechanism for guiding the light from the optical unit to the cartilage and scanning it, and applying a predetermined deformation energy (stress) to the target part. A load mechanism (load device).

Also, a cartilage diagnostic method using the above technique may be constructed. The method includes a step of applying a predetermined deformation energy (load) to the cartilage, a step of displaying the degree of deformation of the cartilage according to the application of the deformation energy (load of load) as a tomographic image by optical coherence tomography, And a step of diagnosing the degree of damage to the cartilage tissue based on the tomographic image.

DESCRIPTION OF THE PREFERRED EMBODIMENTS Embodiments embodying the present invention will be described below in detail with reference to the drawings.
[First embodiment]
FIG. 1 is a diagram schematically illustrating the configuration of the cartilage diagnostic apparatus according to the first embodiment. FIG. 2 is a diagram schematically illustrating a configuration of a diagnostic probe that configures the cartilage diagnostic apparatus. The cartilage diagnostic apparatus of this embodiment can diagnose the degree of cartilage tissue damage (degeneration degree) by applying a predetermined stress to the articular cartilage to be diagnosed and visualizing the degree of cartilage tissue deformation against that stress. It is what. OCT is used for this tomographic visualization.

As shown in FIG. 1, a cartilage diagnostic apparatus 1 includes an optical unit 2 including an optical system using OCT, a diagnostic probe 4 connected to the optical unit 2, and cartilage tissue based on optical interference data obtained by OCT. A control calculation unit 6 for calculating the degree of damage of the cartilage, and a display device 8 for displaying the degree of damage of the cartilage tissue in a manner of visualizing the tomogram. In the illustrated example, an optical system based on a Mach-Zehnder interferometer is shown as the optical unit 2, but a Michelson interferometer or other optical system may be employed.

In this embodiment, SS-OCT (Swept Source OCT) using a wavelength scanning laser as a light source is used as OCT, but TD-OCT (Time Domain OCT), SD-OCT (Spectral Domain OCT) or other OCT is used. You can also. Unlike TD-OCT, SS-OCT does not require mechanical optical delay scanning such as reference mirror scanning, and is preferable in that high time resolution and high position detection accuracy can be obtained.

The optical unit 2 includes a light source 10, an object arm 12, a reference arm 14, and a photodetector 16. Each optical element is connected to each other by an optical fiber. The light emitted from the light source 10 is divided by a coupler 18 (beam splitter), one of which becomes object light guided to the object arm 12 and the other becomes reference light guided to the reference arm 14. The object light guided to the object arm 12 is guided to the probe 4 through the circulator 20, and is irradiated to the cartilage that is the object of diagnosis. This object light is reflected as backscattered light on the surface and cross section of the cartilage, returns to the circulator 20, and is guided to the coupler 22.

On the other hand, the reference light guided to the reference arm 14 is guided to the reflecting mirror 26 via the circulator 24. At this time, the reference light is condensed on the reflecting mirror 26 by the condenser lens 30 through the collimator lens 28. The reference light is reflected by the reflecting mirror 26, returns to the circulator 24, and is guided to the coupler 22. That is, the object light and the reference light are combined (superposed) by the coupler 22, and the interference light is detected by the photodetector 16.

The probe 4 constitutes the object arm 12 of the optical unit 2, and an optical fiber 34 extending from the circulator 20 is inserted into the main body 32 thereof. A syringe needle 36 is attached to the distal end of the main body 32, and the distal end portion of the probe 4 can be inserted into the living body (patient's knee K) through the syringe needle 36. The main body 32 includes an optical mechanism for guiding the light from the optical unit 2 to the cartilage of the knee K to scan, a load mechanism for applying a predetermined load (stress) to the cartilage, and these mechanisms. A drive unit 38 (actuator) for driving is provided. The optical mechanism irradiates light toward the cartilage and guides the reflected light to the object arm 12. Details of each mechanism will be described later.

Interference light combined by the coupler 22 is input to the photodetector 16. The photodetector 16 detects this as an optical interference signal (a signal indicating the intensity of interference light). This optical interference signal is input to the control calculation unit 6 via the A / D converter 40. The A / D converter 40 converts the analog signal output from the photodetector 16 into a digital signal and outputs the digital signal to the control calculation unit 6.

The control calculation unit 6 includes a CPU, a ROM, a RAM, a hard disk, and the like. The hardware and software control the entire optical system of the optical unit 2, drive control of the probe 4, and image output by OCT. For this purpose. The control calculation unit 6 controls the driving of the loading mechanism and the optical mechanism of the probe 4, processes the optical interference signal output from the optical unit 2 based on the driving, and acquires a tomographic image of the cartilage by OCT. Based on the tomographic image data, the degree of damage to the cartilage tissue is calculated by a method described later.

The display device 8 includes, for example, a liquid crystal display, and displays the damage degree of the cartilage tissue calculated by the control calculation unit 6 on the screen in a manner of visualizing the tomography.

As shown in FIG. 2, the probe 4 has a syringe needle 36 as a puncture portion attached to the tip of a cylindrical main body 32. A cylindrical sheath 50 with a bottom is coaxially provided inside the syringe needle 36, and a piezoelectric element 52 (piezo element) capable of driving the sheath 50 in the axial direction is provided inside the main body 32. . The piezoelectric element 52 is supported by a support mechanism (not shown) provided inside the main body 32. The distal end of the sheath 50 is closed, and the rear end opening is fixed to the distal end of the piezoelectric element 52. The distal end surface of the sheath 50 is a contact surface that can contact the surface of the cartilage J of the patient.

The tip of the optical fiber 34 is inserted inside the sheath 50, a GRIN lens 54 is connected to the tip, and a prism mirror 56 is connected to the tip. In other words, the optical fiber 34, the GRIN lens 54, and the prism mirror 56 constitute a light guide integrated in the axial direction. The GRIN lens 54 exhibits a condensing function by continuously changing the refractive index in the lens. The prism mirror 56 is a reflecting surface whose tip surface is inclined with respect to the axis. In this embodiment, the GRIN lens 54 is used as an optical element that adjusts the diameter of the irradiation beam, the focal length, the beam spot diameter, etc., but other optical elements that can perform the same function may be used. . Further, although the prism mirror 56 is employed as an optical element that reflects the beam in a specified direction, other optical elements that can exhibit the same function may be employed.

The light guided to the probe 4 by the optical fiber 34 is emitted from the distal end of the sheath 50 through the GRIN lens 54 and the prism mirror 56, and irradiated to the cartilage J (see the dotted arrow). Thereby, the light reflected on the surface or cross section of the cartilage J is taken into the GRIN lens 54 and guided to the object arm 12 of the optical unit 2 through the optical fiber 34.

The piezoelectric element 52 is provided with an insertion hole 58 penetrating along the axis. The optical fiber 34 is inserted into the insertion hole 58 and introduced into the sheath 50. With such a configuration, a load from the piezoelectric element 52 can be applied to the sheath 50 without applying a mechanical load to the optical fiber 34. By energizing the piezoelectric element 52, the piezoelectric element 52 extends in the axial direction, and the sheath 50 can be driven to the distal end side in the axial direction. This driving force is applied to the cartilage J as a pressing force (compressive stress) by the distal end surface of the sheath 50. That is, the sheath 50 and the piezoelectric element 52 function as a “load device (load mechanism)”.

The sheath 50 is biased rearward by a biasing member (not shown), and is configured to be accommodated inside the syringe needle 36 when the piezoelectric element 52 is not energized. When the piezoelectric element 52 is energized, the distal end surface of the sheath 50 can project from the syringe needle 36 as shown in the figure.

A rotating mechanism 60 (optical rotary joint) for rotating the optical fiber 34 around the axis is provided behind the piezoelectric element 52 inside the main body 32. The rotation mechanism 60 has a rotor fixed coaxially with the optical fiber 34. By energizing the rotation mechanism 60, the optical fiber 34 and the GRIN lens 54 can be rotated around the own axis, and the direction of the prism mirror 56 can be changed. That is, the tip of the optical fiber 34, the GRIN lens 54, the prism mirror 56, and the rotation mechanism 60 function as an “optical mechanism”. By driving the rotation mechanism 60, the object light can be scanned with respect to the cartilage J. Energization control to the piezoelectric element 52 and the rotation mechanism 60 is performed by the control calculation unit 6.

Hereinafter, a calculation processing method until the degree of cartilage tissue damage is presented will be described.
As described above, the object light (reflected light from the cartilage J) that has passed through the object arm 12 and the reference light that has passed through the reference arm 14 are combined and detected by the photodetector 16 as an optical interference signal. The control calculation unit 6 can acquire the optical interference signal as a tomographic image of the cartilage J based on the interference light intensity.

Incidentally, the coherence length l _c which is the resolution in the optical axis direction (depth direction) of OCT is determined by the autocorrelation function of the light source. Here, the coherence length l _c is the half-width of the comprehensive line of the autocorrelation function and can be expressed by the following formula (1).

Here, λ _c is the center wavelength of the laser beam, and Δλ is the full width at half maximum of the laser beam.

On the other hand, the resolution in the direction perpendicular to the optical axis (beam scanning direction) is ½ of the beam spot diameter D based on the light condensing performance of the condensing lens. The beam spot diameter D can be expressed by the following formula (2).

Here, d _B is the beam diameter incident on the condenser lens, f is a focal point of the condenser lens. As described above, the resolution by OCT is limited, but in this embodiment, the diagnosis of cartilage tissue can be performed on a micro scale by introducing a subpixel analysis described later. The details will be described below.

First, a microscale mechanical property detection method using OCT will be described. This detection method calculates a deformation vector distribution by applying a digital cross-correlation method to two OCT tomographic images before and after deformation of a measurement object, and detects a strain rate tensor distribution inside a living tissue on a microscale. It is a technique.

When calculating the deformation vector distribution, the recursive cross-correlation method (Recursive Cross-correlation method), which performs repeated cross-correlation processing, is applied. This is a technique of applying a cross-correlation method by referring to a deformation vector calculated at a low resolution, limiting an exploration area and hierarchically reducing an inspection area. Thereby, a high-resolution deformation vector can be acquired. Further, as a speckle noise reduction method, an adjacent cross-correlation multiplication method (Adjacent cross-correlation Multiplication) that performs multiplication with a correlation value distribution in an adjacent inspection region is used. Then, the maximum correlation value is searched from the correlation value distribution that has been multiplied to increase the SN.

Also, in the micro deformation analysis on the micro scale, the sub-pixel accuracy of the deformation vector is important. For this reason, both the up-stream gradient method using brightness gradient (Up-stream Gradientmethod) and the image deformation method considering image expansion and shear (Image Deformation method) are used together to detect deformation vectors with high accuracy. Is realized. The “windward gradient method” here is a kind of gradient method (optical flow method).

The deformation velocity vector distribution is calculated by time differentiation of the deformation vector thus obtained, and the strain rate tensor distribution is calculated by further spatial differentiation. As will be described later, since the attenuation coefficient of the strain rate tensor of the cartilage tissue can be regarded as equivalent to the damage degree (degeneration degree) of the cartilage, this is visualized as a tomogram as the cartilage damage degree.

Hereinafter, each method will be described in detail. FIG. 3 is a diagram schematically showing a processing procedure by the recursive cross correlation method. FIG. 4 is a diagram schematically showing a processing procedure by subpixel analysis. FIG. 5 is a diagram illustrating a result obtained by each process.

(Recursive cross correlation method)
FIGS. 3A to 3C show processing steps by the recursive cross-correlation method. Each figure shows tomographic images before and after continuous imaging by OCT. The previous tomographic image (Image1) is shown on the left side, and the later tomographic image (Image2) is shown on the right side.

The cross-correlation method is a method for evaluating the local speckle pattern similarity based on the correlation value R _ij based on the following equation (3). Therefore, as shown in FIG. 3A, with respect to the preceding and following OCT images, an inspection area S1 to be inspected with similarity to the previous tomographic image (Image1) is set, and the subsequent tomographic image is displayed. In (Image2), a search area S2 that is a search range of similarity is set.

Here, as the spatial coordinates, the Z axis is set in the optical axis direction, and the X axis is set in the direction perpendicular to the optical axis. f (X _i , Z _j ) and g (X _i , Z _j ) are in the inspection area S1 (N _x × N _z pixels) of the center position (X _i , Z _j ) set in the OCT images before and after the deformation. Represents a speckle pattern.

A correlation value distribution R _{i, j} (ΔX, ΔZ) in the search area S2 (M _x × M _z pixels) is calculated, and pattern matching is performed as shown in FIG. Actually, as shown in the following equation (4), the movement amount U _{i, j} giving the maximum correlation value is determined as the deformation vector before and after the deformation.

Note that f ⁻ and g ⁻ represent average values of f (X _i , Z _j ) and g (X _i , Z _j ) in the inspection region S1.

This method employs a recursive cross-correlation method that increases the spatial resolution by repeating the cross-correlation process while reducing the inspection area S1. In this embodiment, when the resolution is increased, the spatial resolution is doubled. As shown in FIG. 3C, the inspection area S1 is divided into ¼, and the search area S2 is reduced by referring to the deformation vector calculated in the previous hierarchy. Here, the search area S2 is also divided into ¼. By using the recursive cross-correlation method, it is possible to suppress erroneous vectors that occur frequently at high resolution. By performing such recursive cross-correlation processing, the resolution of the deformation vector can be increased as shown in FIG.

In addition to this, by using the following equation (5), a threshold value using the average deviation σ of a total of nine deformation vectors including eight coordinates around the coordinates being calculated is set, and error vectors are removed. Suppresses error propagation associated with recursive processing.

Here, Um represents the median value of the vector quantity, and the coefficient κ serving as a threshold is arbitrarily set.

(Adjacent cross-correlation multiplication)
In this embodiment, an adjacent cross-correlation multiplication method is introduced as a method for determining an accurate maximum correlation value from a highly random correlation value distribution affected by speckle noise. According to the following equation (6), in the adjacent cross-correlation multiplication method, the correlation value distribution Ri, j (Δx, Δz) in the inspection region S1 and Ri + Δi, j (Δx, Δz) with respect to the adjacent inspection region overlapping with the inspection region S1 Multiplication of Ri, j + Δj (Δx, Δz) is performed, and a maximum correlation value is retrieved using a new correlation value distribution R′i, j (Δx, Δz).

This makes it possible to reduce randomness by multiplying correlation values. Since the amount of information of the interference intensity distribution decreases with the reduction of the inspection area S1 described above, it is considered that the appearance of multiple correlation peaks caused by speckle noise causes a deterioration in measurement accuracy. On the other hand, since the amount of movement between adjacent boundaries has a correlation, a strong correlation value remains in the vicinity of the maximum correlation value coordinate. By introducing this adjacent cross-correlation multiplication method, the maximum correlation value peak is clarified, the measurement accuracy is improved, and accurate moving coordinates can be extracted. In addition, by introducing this adjacent cross-correlation multiplication method in each stage of OCT, error propagation is suppressed, and robustness against speckle noise is improved. This makes it possible to calculate the deformation vector and strain distribution with high accuracy even at high spatial resolution.

(Windward gradient method)
4A to 4C show a processing process by subpixel analysis. Each figure shows tomographic images before and after continuous imaging by OCT. The previous tomographic image (Image1) is shown on the left side, and the later tomographic image (Image2) is shown on the right side.

In this embodiment, an upwind gradient method and an image deformation method are used for subpixel analysis. Although the final movement amount is calculated by an image deformation method described later, the windward gradient method is applied prior to the image deformation method because of the problem of convergence of the calculation. The image deformation method and the windward gradient method for detecting the sub-pixel movement amount with high accuracy are applied under the condition of the inspection area size being small and the high spatial resolution. When it is difficult to detect the subpixel movement amount in the image deformation method, the subpixel movement amount is calculated by the windward gradient method.

In sub-pixel analysis, the luminance difference before and after deformation at the point of interest is represented by the luminance gradient and movement amount of each component. For this reason, the sub-pixel movement amount can be determined using the least square method from the luminance gradient data in the inspection region S1. In the present embodiment, when obtaining the brightness gradient, the windward difference method is used which gives the windward brightness gradient before the subpixel deformation.

サ ブ Subpixel analysis strongly depends not only on measurement errors, but also on speckle noise intricately related to tissue scattering effects and polarization characteristics. Furthermore, since the spatial differential coefficient of the movement amount is required for calculating the strain tensor distribution, the subpixel error increases the numerical instability of the differential coefficient calculation, leading to instability of the strain tensor distribution. Various methods exist for the subpixel analysis. In this embodiment, the gradient method for detecting the subpixel movement amount with high accuracy is adopted even under the condition of a small inspection area size and high spatial resolution.

The windward gradient method calculates the movement of the point of interest in the inspection region S1 not only with the pixel accuracy shown in FIG. 4A but also with the sub-pixel accuracy shown in FIG. 4B. Each grid in the figure represents one pixel. Actually, it is considerably smaller than the tomographic image shown in the figure, but for the convenience of explanation, it is shown in large size. This windward gradient method is a method for formulating the change of the luminance distribution before and after the minute deformation by the luminance gradient and the moving amount. When f is the luminance, the following equation (Taylor expansion of the minute deformation f (x + Δx, z + Δz) is performed. 7).

The above formula (7) indicates that the luminance difference before and after the deformation of the attention point is represented by the luminance gradient and the movement amount before the deformation. Since the movement amount (Δx, Δz) cannot be determined only by the above equation (7), the movement amount is considered to be constant in the inspection region S1, and is calculated by applying the least square method.

When calculating the amount of movement using the above equation (7), the luminance difference before and after the movement at each point of interest on the right side can only be obtained uniquely. Therefore, how accurately the luminance gradient is calculated is directly related to the accuracy of the movement amount. In the difference of the brightness gradient, the primary accuracy upwind difference is used. This is because applying high-order differences in differentiating requires a lot of data and is greatly affected when noise is included. In addition, the high-order difference based on each point in the inspection area S1 uses a lot of data outside the inspection area S1, and there is a problem that the amount of movement of the inspection area S1 itself is lost. is there.

When calculating the brightness gradient, it can be considered that the brightness gradient on the windward side before the deformation moves, so that the brightness difference at the point of interest is generated. Therefore, the difference on the windward side is applied before the deformation. The windward here is not the actual movement direction but the direction of the subpixel movement amount with respect to the pixel movement amount, and the upwind is determined by performing parabolic approximation on the maximum correlation value peak. On the contrary, since it can be considered that the luminance difference on the leeward side after the deformation moves in the opposite direction, a difference in luminance at the point of interest occurs, so the difference on the leeward side is applied after the deformation.

2) Two solutions were obtained using the windward difference before deformation and the windward difference after deformation, and the average was taken. Furthermore, when the movement amount does not actually follow the axial direction, the luminance gradient before and after the deformation is not on the same axis as the target point, and it is necessary to obtain a gradient at a shifted position. In other words, when it is desired to obtain the gradient in the X direction, the gradient must be obtained in consideration of movement in the Z direction. Therefore, accuracy is improved by estimating the luminance gradient by luminance interpolation. Basically, a position before deformation (or after deformation) is predicted, and a gradient at the position is obtained by interpolation.

The position of the point of interest before (after) deformation is obtained from the amount of subpixel movement when parabolic approximation is performed. Using the eight coordinates that surround the position of the point of interest, the luminance gradient is calculated from the ratio thereof. Specifically, the following formula (8) is used. The amount of movement was determined by applying the least square method using the luminance gradient thus calculated and the luminance change.

(Image transformation method)
Until the windward gradient method described above, the shape of the inspection region S1 is not changed, and the deformation vector is calculated while keeping the square. However, in reality, it is considered that the inspection region S1 is also deformed in accordance with the deformation of the measurement target. Therefore, it is necessary to introduce an algorithm that takes into account the minute deformation of the inspection region S1 and calculate the deformation vector with high accuracy. There is. For this reason, in this embodiment, an image deformation method is introduced to calculate a deformation vector with sub-pixel accuracy. That is, the cross-correlation is performed between the inspection region S1 before the tissue deformation and the inspection region S1 in consideration of the expansion and contraction and shear deformation after the tissue deformation, and the sub-pixel deformation amount is determined by iterative calculation based on the correlation value. Note that the expansion and contraction and the shear deformation of the inspection region S1 are linearly approximated.

The image deformation method is generally used in a material surface strain measurement method, and is applied to an image obtained by photographing a material surface coated with a random pattern with a high spatial resolution camera. On the other hand, the OCT tomogram not only contains a lot of speckle noise, but especially in a living tissue, the refractive index changes with the flow of the substrate and water in the tissue, so the deformation to the speckle pattern is large. For this reason, a local deformation | transformation generate | occur | produces in a complicated structure | tissue, and when the deformation | transformation of expansion | swelling, contraction, shearing, etc. accompanies inspection area | region S1, the convergence solution of a solution will fall remarkably. Reduction of the examination area S1 in this method is indispensable for detecting local tissue mechanical characteristics. Therefore, in the image deformation method, the amount of deformation obtained by the windward gradient method is adopted as the initial value of the convergence calculation, and further, low robustness is realized even when the inspection area S1 is reduced by bicubic function interpolation of the luminance distribution. ing. In the modification, an interpolation function other than bicubic function interpolation may be used.

More specifically, the calculation is executed according to the following procedure. First, a bicubic function interpolation method is applied to the luminance distribution of the OCT tomogram before tissue deformation, and the luminance distribution is made continuous. The bicubic function interpolation method is a technique for reproducing spatial continuity of luminance information using a convolution function obtained by piecewise approximating a sinc function. Originally, when measuring an image of a continuous luminance distribution, a point spread function depending on the optical system is convolved. Therefore, by performing deconvolution using a sinc function, the original continuous luminance distribution is restored. The When interpolating a discrete uniaxial signal f (x), the convolution function h (x) is expressed by the following equation (9).

Note that the shape of the luminance interpolation function h (x) also needs to be changed depending on the OCT measurement conditions. Therefore, the differential coefficient a at x = 1 of the luminance interpolation function h (x) is variable, and the shape of the luminance interpolation function h (x) can be changed by changing the value of a. In the present embodiment, the value of a is determined based on the verification result by the numerical experiment using the pseudo OCT tomogram. By performing the image interpolation as described above, it is possible to obtain the OCT luminance value at each point in the inspection region S1 in consideration of expansion and contraction and shear deformation.

As shown in FIG. 4C, the inspection region S1 calculated in consideration of expansion and contraction and shear deformation is accompanied by deformation as it moves. Assuming that coordinates (x, z) at integer pixel positions in a certain examination region S1 in the OCT tomogram before tissue deformation move to coordinates (x ^* , z ^* ) after tissue deformation, the values of x ^* , z ^* Is represented by the following formula (10).

Here, Δx and Δz are distances from the center of the inspection region S1 to the coordinates x and z, u and v are deformation amounts in the x and z directions, respectively, and ∂u / ∂x and ∂v / ∂z are x, respectively. The deformation amount in the vertical direction of the inspection region S1 in the z direction, ∂u / ∂z, and ∂v / ∂x are the deformation amounts in the shear direction of the inspection region S1 in the x and z directions, respectively. The Newton-Raphson method is used for the numerical solution, and the correlation value derivative with 6 variables (u, v, ∂u / ∂x, ∂u / ∂z, ∂v / ∂x, ∂v / ∂z) is 0. That is, iterative calculation is performed so as to obtain the maximum correlation value. In order to improve the convergence of the iterative calculation, the sub-pixel movement amount obtained by the windward gradient method is used as the initial movement amount in the x and z directions. When the Hessian matrix for the correlation value R is H and the Jacobian vector for the correlation value is ▽ R, the update amount ΔPi obtained in one iteration is expressed by the following equation (11).

For convergence judgment, it is used that the asymptotic solution obtained at any time by iterative calculation is sufficiently small in the vicinity of the convergence solution. However, in a region where the speckle pattern changes drastically, a correct convergence solution may not be obtained because it cannot be followed by linear deformation. In this case, the present embodiment employs the sub-pixel movement amount obtained by the windward gradient method.

A deformation velocity vector distribution as shown in FIG. 5B can be calculated by differentiating the deformation vector of subpixel accuracy obtained as described above with respect to time. A strain rate tensor can be calculated by spatially differentiating the deformation rate vector distribution.

(Time-space moving least squares method)
For the calculation of the strain rate tensor, a moving least square method (hereinafter referred to as “MLSM”) is used. MLSM is a technique that enables smooth calculation of a differential coefficient while smoothing a movement amount distribution. The square error equation used in MLSM is expressed by the following equation (12).

In the above equation (12), parameters a to k that minimize S (x, z, t) are obtained. That is, the following equation (13) is adopted as a three-variable quadratic polynomial in the horizontal direction x, the depth direction z, and the time direction t as an approximation function. Based on the least square approximation, an optimal derivative is calculated from the following equation (14) and smoothed.

By using this derivative, the strain rate tensor shown in the following formula (15) can be calculated. fx and fz indicate the strain increment of each axis, and the strain rate is calculated from the amount of change over time.

Next, the cartilage diagnosis process performed using each of the above-described methods will be described. 6 and 7 are diagrams illustrating a processing procedure of the cartilage diagnosis processing.
In this embodiment, the above-described methods are applied to the cartilage diagnostic apparatus 1. As shown in FIG. 6, the control calculation unit 6 sequentially stores the deformation velocity vector distribution calculated as described above (FIG. 6A), and is obtained by spatially differentiating each deformation velocity vector distribution. The strain rate tensor distribution is sequentially stored (FIG. 6B). As shown in FIG. 7A, the strain rate tensor distribution can be displayed in the form of a tomographic image so as to be identifiable for each predetermined range. In the figure, the region where the strain rate is negative indicates the compression region, and the positive region indicates the expansion region.

Then, when the temporal change in strain rate for each coordinate in the tomographic image is calculated, a calculation result as shown in FIG. 7B is obtained. In the figure, the horizontal axis indicates the passage of time, and the vertical axis indicates the strain rate. When smoothing processing is performed on the time change of the strain rate calculated in this way, and the strain rate attenuation coefficient is calculated for each coordinate, an attenuation coefficient distribution as shown in FIG. 7C can be obtained. This tomographic distribution of the attenuation coefficient can be evaluated as corresponding to the distribution of the degree of damage of the cartilage tissue, as can be seen from the experimental results described later. Therefore, the control calculation unit 6 causes the display device 8 to visualize the tomographic distribution by making the tomographic distribution of the attenuation coefficient equivalent to the distribution of the degree of damage of the cartilage tissue.

In this embodiment, the fault distribution of the attenuation coefficient itself is presented as the degree of cartilage tissue damage (degeneration degree). Therefore, a region with a large attenuation coefficient is indicated as a region with a high degree of damage (modification degree), and a region with a low attenuation coefficient is indicated as a region with a low degree of damage (modification degree). In addition, the hatching part in a figure shows the area | region (area | region with low precision) which has a small strain rate and shakes to a positive / negative value, and can be disregarded in cartilage diagnosis.

Next, the results of experiments conducted to verify that the fault distribution of the attenuation coefficient corresponds to the distribution of the degree of damage of the cartilage tissue will be described. FIG. 8 is a diagram schematically illustrating the configuration of the experimental apparatus. FIG. 9 is a diagram illustrating a test piece used in the experiment. FIG. 10 is a graph showing a state of a load applied to the test piece during the experiment. In FIG. 8, the same components as those shown in FIG.

As shown in FIG. 8, this experimental apparatus is configured by connecting a compression tester 70 to the optical unit 2 shown in FIG. The compression tester 70 applies a compressive load to the test piece W instead of the probe 4 shown in FIG. As the test piece W, cartilage of the knee joint of a pig is used. The compression tester 70 includes a glass window 72 that contacts the test piece W, a metal indenter 74 that sandwiches the test piece W between the glass window 72, a linear actuator 76 that drives the metal indenter 74 in the axial direction, and the test piece W. The load cell 78 etc. which detect the load loaded on is provided.

In this experiment, the control arithmetic unit 6 is composed of two personal computers 80 and 82, one of which drives the optical unit 2 and OCT arithmetic processing, and the other drives the compression tester 70. A test probe 84 is connected to the object arm 12. The test probe 84 includes an optical mechanism capable of scanning the object light toward the test piece W instead of the probe 4 shown in FIG.

9A shows an installation mode of the test piece W to the compression tester 70, FIG. 9B shows a side view and the size of the test piece W, and FIG. The front view and its size are shown. As a test piece W, a cartilage tissue sample of a femoral condyle of a knee joint of a pig about 6 months old was collected in a cylindrical shape having a diameter of 5.5 mm. In addition, in order to diagnose and evaluate the mechanical properties of degenerated cartilage, simulated degenerated cartilage was prepared by enzyme treatment and compared with normal cartilage. For the enzyme treatment, collagenase type 2 (C6885, Sigma Ardrich) that degrades collagen fibers maintaining the elasticity of articular cartilage was used. Normal articular cartilage is placed in 30 [unit / ml] collagenase solution prepared by dissolving collagenase in phosphate buffered saline (10010, Invitrogen) at pH 7.4, and 0.5 [37] at 37 [° C.]. The treatment is performed with a time setting of 1 hour, 2 hours, and collagen fiber degeneration is simulated as the degree of OA degeneration.

Returning to FIG. 8, in this experiment, the lower bone side of the test piece W was fixed to the metal indenter 74 and pressed against the glass window 72 installed in the compression tester 70 to cause compression deformation in the cartilage tissue. An OCT tomographic image during the stress relaxation test was continuously acquired at 27.03 [flame / sec] by implanting an OCT beam through the glass window 72. In this compression experiment, the compression amount was fixed at 10 [%] strain, and the compression speed was fixed at 0.1 [% / sec]. The compression load loading time was 100 [sec], and after reaching 10 [%] strain, the behavior of the cartilage tissue during stress relaxation was observed. In addition, the recursive cross-correlation method is set to two layers for reducing the inspection area from 9 [pixel] to 5 [pixel], and the subpixel accuracy is set at a vector density of 6.5 × 15.62 [μm]. The deformation velocity vector distribution of was calculated.

FIG. 10 shows stress time series data detected by the load cell 78 in this stress relaxation test. In the figure, the horizontal axis indicates time [sec] from the start of load application, and the vertical axis indicates compressive stress [MPa]. As shown in the figure, the stress increases during compression (0 to 100 seconds), and stress relaxation occurs after the end of compression (100 seconds). It was observed that an equilibrium state was reached around 1200 seconds. In addition, since the stress is attenuated exponentially in this way, it is understood that the cartilage is a viscoelastic body whose relaxation elastic modulus changes with time.

FIG. 11 and FIG. 12 are diagrams in which a change in strain rate distribution over time by a stress relaxation test is visualized as a fault. FIG. 11 shows the calculation results for normal cartilage, (A) is 0.31 [sec] from the start of loading, (B) is 50.35 [sec], (C) is 100.29 [sec], ( D) is 101.96 [sec], (E) is 110.95 [sec], (F) is 120.12 [sec], (G) is 130.52 [sec], and (H) is 140.51. [sec] and (I) show the result of 150.49 [sec].

On the other hand, FIG. 12 shows the calculation results for the simulated degenerated cartilage, (A) is 1.07 [sec] from the start of loading, (B) is 50.07 [sec], and (C) is 100.99 [sec]. ], (D) is 101.99 [sec], (E) is 110.98 [sec], (F) is 120.97 [sec], (G) is 130.96 [sec], (H) is 140.95 [sec], (I) shows the result of 150.94 [sec]. That is, (A) and (B) in each figure indicate that compression is in progress, and (C) to (I) indicate that stress is being relaxed (see FIG. 10).

11 and 12, it can be confirmed that a large compressive strain rate is generated in the surface layer of the test piece W. This region corresponds to the surface layer region (15% from the surface) of the cartilage, and corresponds to the region where the deformation rate vector changes rapidly. On the other hand, 15 to 70% from the cartilage surface is called an intermediate layer region and corresponds to a region where the detected strain rate is small and a rapid change in the deformation rate vector is not observed. It is known that the orientation characteristics of collagen fibers are greatly different between the surface layer and the intermediate layer. In particular, since the surface layer is horizontally oriented on the cartilage surface, the rigidity in the compression direction is low, and it is considered that the concentration of the compressive strain rate occurred.

Further, when FIG. 11 and FIG. 12 are compared, the remaining strain rate can be confirmed even at 150.49 [sec] in normal cartilage, whereas the strain rate is almost observed at 130.96 [sec] in simulated degenerated cartilage. You can see that there is no. Therefore, it can be seen that the degree of degeneration of articular cartilage can be evaluated by focusing on the relaxation time and comparing. From this, it is considered that the degree of cartilage degeneration can be evaluated by comparing the attenuation coefficient calculated from the change in strain rate with time. In this embodiment, paying attention to this point, as shown in FIG. 7, the distribution of the attenuation coefficient is calculated from the strain rate tensor distribution. That is, the tomographic distribution of the strain rate attenuation coefficient is visualized on the screen as the equivalent of the distribution of the damage degree of the cartilage tissue, thereby enabling cartilage diagnosis.

FIG. 13 is a view showing a comparison of attenuation coefficients of normal cartilage and simulated degenerated cartilage. The horizontal axis in the figure indicates the enzyme treatment time (H: hour), and the vertical axis indicates the attenuation coefficient. The results of 0H, 0.5H, 1H, and 2H are shown as the enzyme treatment time. Note that 0H indicates zero enzyme treatment time, that is, normal cartilage. The result shown is a spatial average of the attenuation coefficient distribution in the surface layer region, and a logarithmic process is applied to the relaxation behavior from the start of stress relaxation to 50 seconds.

From the results shown in the figure, it can be seen that the attenuation coefficient increases as the enzyme treatment time increases. In collagenase enzyme treatment, with the degradation of collagen fibers that govern elasticity, the proteoglycans that govern viscosity also occur. For this reason, in addition to the increase in porosity due to collagen fiber degradation, the water retention property is significantly increased due to the decrease in water retention property due to proteoglycan detachment, and the viscosity property is considered to be impaired. Since the increase in water permeability means failure of the load support mechanism due to viscosity, it is considered that the damping coefficient of the relaxation elastic coefficient decreases as shown in the figure. As shown in the figure, there was a significant difference in the results of any enzyme treatment time with respect to normal cartilage, and it was possible to detect degeneration. The simulated degenerated cartilage sample used in this experiment has no change in appearance as compared with normal cartilage, and is denatured to the extent that it cannot be judged by visual diagnosis with X-rays or CT. Therefore, it can be seen that the present method can diagnose early knee osteoarthritis that could not be diagnosed.

Next, the flow of specific processing executed by the control calculation unit 6 will be described.
FIG. 14 is a flowchart showing a flow of cartilage diagnosis processing executed by the control calculation unit 6. This cartilage diagnosis process is started in a state where the tip of the probe 4 is inserted into the patient's knee joint by a doctor or the like. The control calculation unit 6 drives the drive unit 38 to execute the stress load and stress relaxation processing shown in FIG. 10, while executing the processing shown in FIG.

That is, the control calculation unit 6 reads two consecutive OCT tomographic images I (x, z, t) and I (x, z, t + Δt) taken at different times (S10). Subsequently, processing by a recursive cross correlation method is executed. Here, first, a cross-correlation process is executed at the minimum resolution (inspection area of the maximum size) to obtain a correlation coefficient distribution (S12). Subsequently, a product of adjacent correlation coefficient distributions is calculated by the adjacent cross correlation multiplication method (S14). At this time, error vectors are removed by a spatial filter such as a standard deviation filter (S16), and interpolation of the removal vectors is executed by a least square method or the like (S18). Subsequently, the cross correlation process is continued by increasing the resolution by reducing the inspection area (S20). That is, the cross-correlation process is executed based on the reference vector at the low resolution. If the resolution at this time is not the predetermined maximum resolution (N in S22), the process returns to S14.

Then, the processing of S14 to S20 is repeated, and when the cross-correlation processing at the highest resolution is completed (Y in S22), sub-pixel analysis is executed. That is, based on the distribution of the deformation vector at the maximum resolution (the inspection area of the minimum size), the subpixel movement amount by the windward gradient method is calculated (S24). Based on the sub-pixel movement amount calculated at this time, the sub-pixel deformation amount by the image deformation method is calculated (S26). Subsequently, the error vector is removed by filtering using the maximum cross-correlation value (S28), and interpolation of the removal vector is executed by the least square method or the like (S30). Then, time differentiation of the deformation vector obtained in this way is performed, and the tomographic distribution U (x, z, t), W (x, z, t) of the deformation velocity vector is calculated for the tomographic image, and the process goes to S10. Return. Here, U (x, z, t) is the z-direction component of the deformation velocity vector distribution, and W (x, z, t) is the x-direction component of the deformation velocity vector distribution. The above processing is also executed for subsequent tomographic images (N in S34). Then, when the tomographic distribution of the deformation vector is obtained for the set number of tomographic images (Y in S34), a modification degree calculation presenting process for cartilage diagnosis is executed (S36).

FIG. 15 is a flowchart showing in detail the modification degree calculation presenting process of S36 in FIG. In the modification degree calculation presenting process, the control calculation unit 6 executes the smoothing of the deformation velocity vector distribution calculated in S34 by the spatiotemporal movement least square method (S40). Then, a strain rate tensor is calculated by performing spatial differentiation on the smoothed deformation rate vector (S42). Then, logarithmic transformation is performed on the strain rate time-series data (S44), and the strain rate attenuation coefficient is calculated by the least square method (S46). The distribution of the attenuation coefficient obtained in this way is regarded as equivalent to the degree of cartilage tissue damage, and a tomographic visualization is made on the display device 8 (S48). That is, for example, the degree of cartilage degeneration is presented in the manner shown in FIG.

As described above, according to the present embodiment, a diagnosis target is obtained by an arithmetic process using the correspondence between the strain rate change (attenuation coefficient) due to stress relaxation and the degree of cartilage tissue damage. The degree of cartilage damage is visualized as a tomogram. For this reason, it becomes possible for a doctor or the like to easily perform cartilage diagnosis by confirming the image visualized by the tomography. That is, cartilage diagnosis using OCT can be more practically used.

[Second Embodiment]
Next, a second embodiment of the present invention will be described. The present embodiment is almost the same as the first embodiment except that a dynamic viscoelasticity method is applied instead of the stress relaxation method as a stress loading method. For this reason, the same components as those in the first embodiment are denoted by the same reference numerals as necessary, and the description thereof is omitted as appropriate.

In this embodiment, a predetermined vibration load (dynamic force, dynamic load) is generated by the probe 4 shown in FIGS. That is, when the probe 4 is energized, the piezoelectric element 52 vibrates at a predetermined excitation frequency, and a vibration load is applied to the cartilage J. Thereby, the dynamic viscoelasticity of the cartilage J is evaluated, and the degree of damage (degree of degeneration) is visualized as a tomogram.

FIG. 16 is a flowchart showing in detail the modification degree calculation presenting process according to the second embodiment. This process is executed in place of the process shown in FIG. 15 of the first embodiment. The control calculation unit 6 drives the drive unit 38 to apply a dynamic load (dynamic force, vibration load) to the cartilage J, while executing the processing shown in FIGS. 14 and 16.

In this embodiment, in the process of S24 shown in FIG. 14, the luminance distribution is formulated using the following equation (16) instead of the above equation (7).
δf / δt = δf / δx · Δx + δf / δy · Δy + δf / δz · Δz (16)

Then, as shown in FIG. 16, in the modification degree calculation presenting process, the control calculation unit 6 executes the smoothing of the deformation velocity vector distribution calculated in S34 by the spatio-temporal moving least square method (S240). Then, Fourier transformation is performed on the smoothed deformation velocity vector (S242). As a result, a bandpass filter process for extracting only the component synchronized with the excitation frequency ω of the probe 4 in the deformation velocity vector is executed (S244). Thereafter, the fault distribution of the deformation velocity vectors U _ω (x, z, t) and W _ω (x, z, t) of the excitation frequency component is calculated by inverse Fourier transform (S246). Here, U _ω (x, z, t) is the z-direction component of the deformation speed vector, and W _ω (x, z, t) is the x-direction component of the deformation speed vector. W _ω (x, z, t) can be expressed, for example, in the form of the following formula (17).
_{W ω (x, z, t} ) = A ω (x, z) sin (ωt + Φ ω (x, z)) ··· (17)

Subsequently, the tomographic distribution of the normalized amplitude A ′ _ω (x, z) and the phase difference Φ _ω (x, z) of the deformation velocity vector related to the excitation frequency component is calculated from the displacement information of the piezoelectric element 52 (S248). ). Here, the normalized amplitude A ′ _ω (x, z) is obtained by dividing the amplitude A _ω (x, z) by the amplitude A _pzt of the displacement information W _pzt (t) = A _pzt sin (ωt) of the piezoelectric element 52. Can be obtained. In the modification, the normalized amplitude A ′ _ω (x, z) and the phase difference Φ _ω (x, z) may be directly calculated by the least square method.

Subsequently, the dynamic viscoelasticity of the cartilage is calculated based on at least one tomographic distribution of the normalized amplitude A ′ _ω (x, z) and the phase difference Φ _ω (x, z) (S250). Then, based on the calculated dynamic viscoelasticity, the damage degree of the cartilage tissue is visualized as a tomogram (S48).

Various methods for evaluating the dynamic viscoelasticity can be considered. When the viscoelasticity evaluation is performed based on the fault distribution of the normalized amplitude A ′ _ω (x, z), it can be evaluated as follows. That is, the dynamic viscoelastic modulus becomes small at a location where the normalized amplitude A ′ _ω (x, z) is large. That is, the normalized amplitude A ′ _ω (x, z) is relatively large in the surface layer where the rigidity in the depth direction of the cartilage is low, and relatively small in the middle layer. As the cartilage degeneration progresses, the rigidity in the depth direction of the surface layer is further reduced, and thus the normalized amplitude A ′ _ω (x, z) is considered to be further increased.

When viscoelasticity evaluation is performed based on the fault distribution of the phase difference Φ _ω (x, z), phase connection processing (unwrapping processing) may be performed. In this case, it can be evaluated as follows. That is, the phase difference Φ _ω (x, z) in the depth direction is zero on the surface of the cartilage, and gradually increases on the surface layer toward the deep part due to the viscoelastic effect of the surface layer. That is, the phase is delayed. In denatured cartilage, the viscosity characteristics of the surface layer are impaired, so that the phase difference Φ _ω (x, z) in the depth direction becomes small (that is, it exhibits elastic behavior).

The degree of degeneration can also be diagnosed by the magnitude of the spatial differential (∂Φ _ω / ∂z) of the phase difference Φ _ω (x, z). A portion where ∂Φ _ω / ∂z is large has a low degree of modification (that is, normal), and a portion where ∂Φ _ω / ∂z is small is said to have a high degree of modification (that is, the viscosity effect is impaired). Based on such knowledge, it is possible to visualize the degree of damage to the cartilage tissue.

The preferred embodiments of the present invention have been described above, but the present invention is not limited to the specific embodiments, and various modifications can be made within the scope of the technical idea of the present invention. Nor.

In the first embodiment, an example is shown in which the attenuation coefficient of the strain rate is calculated and the tomographic visualization is equivalent to the degree of cartilage damage. In the modified example, the deformation rate attenuation coefficient may be calculated, and the tomographic visualization may be made equivalent to the degree of cartilage damage.

In the second embodiment, the example in which the viscoelasticity is evaluated from the normalized amplitude A ′ _ω (x, z) and the phase difference Φ _ω (x, z) has been shown. However, the normalized amplitude A ′ _ω (x, z A value expressing a viscoelastic physical quantity may be newly defined from z) and the phase difference Φ _ω (x, z), and viscoelasticity evaluation may be performed based on the value.

In the second embodiment, the viscoelasticity of the cartilage is evaluated based on the tomographic distribution of the deformation velocity vector in the modification degree calculation presenting process shown in FIG. In the modified example, the tomographic distribution of the strain rate tensor may be calculated, and the viscoelasticity of the cartilage may be similarly evaluated based on the calculated tomographic distribution. This strain rate tensor can be calculated by spatially differentiating the deformation rate vector.

Although not clearly shown in the second embodiment, a sensor such as a load cell may be attached to the probe 4 to measure the load applied to the cartilage. Thereby, the complex elastic modulus (storage elastic modulus, loss elastic modulus, complex viscosity, etc.) of the cartilage tissue can be calculated, and the normalized amplitude A ′ _ω (x, z) and the phase difference Φ _ω (x, z ) Makes it easy to define physical quantities for evaluating viscoelasticity.

Although not described in the above embodiment, an amplitude value of strain rate and a time delay (phase delay) of strain rate may be adopted as the “mechanical feature amount”. In other words, taking the dynamic viscoelastic method described above as an example, the strain rate fluctuates between a positive value and a negative value in the process of repeatedly applying and relaxing stress. The strain rate also varies so as to follow the variation of the load. In this regard, the magnitude (amplitude value) of the strain rate fluctuation and the followability (responsiveness) of the strain rate fluctuation to the load fluctuation tend to change corresponding to the degree of cartilage damage. Specifically, as the degree of cartilage damage progresses, the amplitude value of the strain rate increases and the time delay (phase delay) tends to decrease. Therefore, the fault distribution may be calculated with respect to the amplitude value of the strain rate and the time delay (phase delay). The fault distribution with the amplitude or time delay (phase delay) may be visualized as equivalent to the distribution of the damage degree of the cartilage tissue.

Alternatively, the median value of the periodic fluctuation of the strain rate may be adopted as the “mechanical feature value”. In other words, taking the dynamic viscoelastic method described above as an example, the strain rate fluctuates between a positive value and a negative value in the process of repeatedly applying and relaxing stress. The center of the fluctuation of the strain rate tends to deviate slightly from zero due to the balance between the viscous force and the elastic force. And the amount of deviation tends to change corresponding to the degree of cartilage damage. Therefore, the fault distribution may be calculated for the amount of deviation from zero of the variation center (median value) of the strain rate. The tomographic distribution of the deviation amount may be visualized as a tomogram equivalent to the distribution of the degree of damage of the cartilage tissue.

In the above embodiment, a load mechanism that applies stress by contact with a piezoelectric element or the like is illustrated as a load device for applying predetermined deformation energy (load) to cartilage. In the modified example, a load device that applies a load (excitation force) to the cartilage in a non-contact manner using ultrasonic waves (sound pressure), photoacoustic waves, electromagnetic waves, or the like may be employed.

In addition, this invention is not limited to the said Example and modification, A component can be deform | transformed and embodied in the range which does not deviate from the summary. Various inventions may be formed by appropriately combining a plurality of constituent elements disclosed in the above-described embodiments and modifications. Moreover, you may delete some components from all the components shown by the said Example and modification.

Claims (8)

1. A cartilage diagnostic apparatus for diagnosing articular cartilage,
   An optical unit including an optical system using optical coherence tomography;
   An optical mechanism that is connected to the optical unit and has a tip that can be inserted into a joint space, guides the light from the optical unit to the cartilage, and scans the cartilage, and applies predetermined deformation energy to the cartilage. A probe including a load device for
   Controls the driving of the load device and the optical mechanism, processes the optical interference signal output from the optical unit based on the driving, and has a predetermined mechanical characteristic amount associated with deformation inside the cartilage by applying deformation energy A control calculation unit that calculates the change in association with the tomographic position of the cartilage and calculates the degree of damage to the cartilage tissue based on the change in the mechanical feature amount;
   A display device for displaying the degree of damage of the cartilage tissue in a manner to visualize a tomography;
   A cartilage diagnostic apparatus comprising:
2. The control calculation unit is
   Relieving stress after applying a predetermined load to the cartilage by driving the load device,
   Based on tomographic image data taken by the optical unit before and after the stress relaxation, a deformation velocity vector corresponding to a cartilage tomographic position is calculated as the mechanical feature amount, and an attenuation coefficient obtained from the change of the deformation velocity vector The cartilage diagnostic apparatus according to claim 1, wherein the tomographic distribution is visualized as a tomogram equivalent to the distribution of the degree of damage of the cartilage tissue.
3. The control calculation unit calculates a strain rate tensor by spatially differentiating the deformation rate vector calculated as the mechanical feature value, and assumes that the fault distribution of the strain rate attenuation coefficient is equivalent to the distribution of the damage degree of the cartilage tissue. The cartilage diagnostic apparatus according to claim 2, wherein a tomographic visualization is performed.
4. The said control calculating part calculates the said deformation | transformation speed vector with sub pixel precision based on the cross correlation of the tomographic image data acquired sequentially, the displacement of a brightness | luminance gradient, and the deformation | transformation of a brightness pattern. The cartilage diagnostic apparatus according to the description.
5. The control calculation unit is
   A predetermined dynamic force is applied to the cartilage by driving the load device;
   Based on the tomographic image data taken by the optical unit during the vibration of the cartilage due to the dynamic force, a deformation velocity vector at each tomographic position of the cartilage is calculated as the mechanical feature amount, and the deformation velocity vector is added. Calculates the viscoelastic tomographic distribution of cartilage based on at least one of the amplitude and phase difference of the specific deformation velocity vector, which is the vibration frequency component, and visualizes the tomogram as the viscoelastic tomographic distribution is equivalent to the distribution of the degree of damage of the cartilage tissue. The cartilage diagnostic apparatus according to claim 1, wherein:
6. The control calculation unit is
   A predetermined dynamic force is applied to the cartilage by driving the load device;
   Based on the tomographic image data taken by the optical unit during the vibration of the cartilage by the dynamic force, the strain rate tensor at each tomographic position of the cartilage is calculated as the mechanical feature amount, and the strain rate tensor is added. Calculates the viscoelastic tomographic distribution of cartilage based on at least one of the amplitude and phase difference of the specific strain rate tensor, which is the vibration frequency component, and visualizes the tomographic fault as equivalent to the distribution of the degree of damage of the cartilage tissue. The cartilage diagnostic apparatus according to claim 1, wherein:
7. The load device includes a sensor capable of measuring a force or load applied to cartilage,
   The control calculation unit is
   A viscoelastic tomographic distribution of cartilage is calculated based on at least one of the amplitude and phase difference of the excitation frequency component of the calculated mechanical feature quantity and a measured value by the sensor, and the viscoelastic tomographic distribution is calculated as the cartilage tissue. 7. The cartilage diagnostic apparatus according to claim 5 or 6, wherein a tomographic visualization is made equivalent to the distribution of the degree of damage.
8. A diagnostic probe that enables diagnosis of a target site by being connected to an optical unit including an optical system that uses optical coherence tomography,
   A tip portion configured to be able to contact the target portion;
   An optical mechanism for guiding the light from the optical unit to the cartilage and scanning;
   A load device for applying predetermined deformation energy to the target part;
   A diagnostic probe comprising:

Images