JP2020096763A - Image processing apparatus and magnetic resonance imaging apparatus - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a new analysis method different from a diffusion distribution analysis method using a model or approximation. (EN) A method capable of presenting the spatial distribution of diffusion coefficients in an easy-to-understand format with high tissue separation ability for complex tissues. In an image having spatial information such as a diffusion coefficient image, a spatial distribution of a physical quantity that is a pixel value of the image is analyzed. At that time, only the spatial information inherent in the pixel value is extracted, and the mutual relationship (spatial relationship) between the plurality of spatial information is grasped. Then, using this spatial relationship, a texture analysis method is applied to each voxel to calculate a texture index, and an image having this index as a voxel value is created. The image created in this way is an image that briefly shows the distribution of physical quantities, for example, the degree of uniformity, the degree of disorder, and the local correlation. [Selection diagram] Figure 2

Description

The present invention relates to a technique for processing an image acquired by a medical imaging device, and particularly to a technique for processing an image acquired by a magnetic resonance imaging device.

A magnetic resonance imaging (hereinafter abbreviated as MRI) apparatus applies a radio frequency magnetic field (RF magnetic field) and a gradient magnetic field to a subject placed in a static magnetic field, thereby obtaining a nuclear magnetic resonance signal generated from the subject. This is an apparatus for imaging and obtains an image emphasizing the physical quantity of the subject tissue or blood flow that characterizes the nuclear magnetic resonance signal by changing the conditions of application of the RF pulse and the gradient magnetic field pulse, and further It can be derived. Such physical quantity is important information as a clue to know degeneration of tissues and lesions in blood vessels.

For example, one of the imaging methods using MRI is a diffusion weighted image (DWI: diffusion weighted image). The diffusion weighted image is an image that emphasizes the self-diffusion of water molecules contained in the biological tissue. It is known that lesions can be visualized immediately after the onset of acute cerebral infarction, and a contrast different from that of T1-weighted images and T2-weighted images is shown. Generally, an imaging method is used in which a signal intensity of randomly moving nuclear spins is lowered by applying a strong gradient magnetic field called MPG (motion probing gradient) pulse. Since the decrease in signal intensity is caused by the nuclear spins that diffuse in the application direction of MPG, it is possible to obtain diffusion information in any direction by controlling the application direction of MPG. Further, the diffusion emphasis degree can be adjusted by a diffusion factor (b value) which is a parameter relating to the application intensity of MPG and time, and the higher the b value, the higher the degree of diffusion emphasis image can be acquired. Furthermore, the diffusion-weighted image can be analyzed to calculate an index relating to diffusion such as a diffusion coefficient. The diffusion coefficient calculated for each application axis is called an apparent diffusion coefficient (ADC: Appearance Diffusion Coefficient). Further, in DTI (diffusion tensor imaging) assuming a three-dimensional elliptical diffusion model of normal distribution, it is possible to calculate a mean diffusion coefficient (MD) and anisotropy ratio (FA).

The MD is an index for knowing the average size of diffusion, and there is no information on the direction of diffusion. On the other hand, FA contains information on diffusion distribution and is widely used as a method for analyzing the structure of the nerve track of white matter, which is a highly anisotropic tissue structure. On the other hand, the FA value decreases in the white matter region where the nerve fibers have uneven kissing or crossing structures or gray matter with small diffusion anisotropy, so in the FA that has only simple diffusion distribution information, It becomes difficult to grasp the structure.

To solve this problem, a plurality of fibers are obtained by a method (for example, Non-Patent Document 1) of performing Funk-Radon transformation on data acquired by applying MPG in multiple directions to obtain an orientation distribution function (ODF). Analytical methods have been proposed that can be applied to a structure in which the two intersect. In addition, in order to estimate specific information such as density and direction variation of axons and dendrites in the brain, a diffusion model consisting of three compartments of cerebrospinal fluid component, intracellular diffusion and extracellular diffusion is assumed. A method of performing the analysis (for example, Non-Patent Document 2) has been proposed.

On the other hand, various methods have been proposed in which an image captured by a medical imaging apparatus is processed by using an image processing technique to create an image and information useful for diagnosis. For example, Patent Document 1 discloses a technique of dividing a region estimated to be a lesion in a medical image and extracting a feature amount of each divided region. To extract the feature amount, a histogram or a texture analysis result is calculated as a feature amount vector. Texture analysis is a method that digitizes spatial patterns of images and performs image classification and segmentation.

JP, 2016-7270, A

Maxime Descoteaux et al., "Regularized, Fast, and Robust Analytical Q-Ball Imaging", Magnetic Resonance in Medicine 58:497-510 (2007). Hui Zhang et al., "NODDI: Practical in vivo neural orientation dispersion and density imaging of the human brain", NeuroImage 61 (2012) 1000-1016.

Since the information provided by the method disclosed in Non-Patent Document 1 is given a complicated ODF shape for each voxel, it is difficult to grasp the degree of diffusion in that region at a glance. Further, the method disclosed in Non-Patent Document 2 can grasp the states of axons and dendrites, but it is difficult to interpret the calculated values for the region that does not fit the set model, and the image It is not possible to analyze diffusion information with a common interpretation throughout.

On the other hand, as disclosed in Patent Document 1, the texture analysis is a method generally applied to a two-dimensional image, and can be applied to an image having a diffusion index as a value such as an MD image or an FA image. .. This makes it possible to clarify the difference in tissues characterized by MD and FA, and to provide information useful for segmentation of tissues and discrimination between normal tissues and lesion sites. However, since the analysis combined with the DTI for calculating the diffusion index from the measurement data and Non-Patent Documents 1 and 2 is required, the problem of each diffusion analysis method is still unsolved.

It is an object of the present invention to provide a technique capable of presenting the characteristics of the distribution of physical quantities as voxel values of an image in an easy-to-understand manner when imaging a physical quantity having a spatial distribution within a voxel. In particular, it is an object to provide a technique capable of presenting the distribution of physical quantities such as diffusion coefficient in a form that is easy to diagnose.

In order to solve the above-mentioned problems, the present invention extracts only spatial information that is inherent in pixel values and grasps the mutual relationship (spatial relationship) of a plurality of spatial information. Then, using this spatial relationship, a texture analysis method is applied to each voxel to calculate a texture index, and an image having this index as a voxel value is created. The image created in this way is an image that briefly shows the distribution of physical quantities in voxels, for example, the uniformity, the degree of randomness, the local correlation, etc., and is an image useful as a diagnostic assistance image.

That is, the image processing apparatus of the present invention includes a reception unit that receives a plurality of image data having a physical quantity as a pixel value and the physical quantity has spatial information, and a mutual position that extracts the relationship between the spatial information in each of the plurality of image data. A calculation unit and a texture analysis unit that analyzes the distribution of the physical quantity and calculates a texture index using the spatial relationship extracted by the mutual position calculation unit. The image processing apparatus of the present invention may further include an image generation unit that generates an image in which the texture index calculated for each voxel by the texture analysis unit is used as a pixel value.

Further, the MRI apparatus of the present invention includes a measuring unit that measures a nuclear magnetic resonance signal using a plurality of MPGs having different application axes, and an image reconstructing unit that creates an image using the nuclear magnetic resonance signal measured by the measuring unit. And an image processing unit that processes the image created by the image reconstructing unit, and the image processing unit has the functions of the image processing apparatus.

According to the present invention, the spatial distribution of the physical quantity included in the data can be presented as a texture index and further as an image representing the texture index without using a model or approximation. As a result, the spatial distribution of the physical quantity can be easily grasped, the classification function of the target tissue is improved, and useful information for diagnosis can be presented.

FIG. 1 is a block diagram showing the overall configuration of an image processing device according to a first embodiment. (A) is a figure which shows the outline of a process of the image processing apparatus of 1st embodiment, (B) is a figure which shows the conventional process. The block diagram which shows the whole structure of the MRI apparatus of 2nd embodiment. The functional block diagram of the computer of the MRI apparatus of 2nd embodiment. The flowchart which shows operation|movement of the MRI apparatus of 2nd embodiment. The figure which shows an example of the pulse sequence of diffusion weighted imaging. (A)-(C) is a figure explaining a triangle surface mesh as an example of spatial relation. The figure explaining GLCM used by general texture analysis. The figure which shows the process sequence of the texture analysis of 2nd embodiment. The figure explaining the processing procedure of the texture analysis of 2nd embodiment. The flowchart which shows the detail of the processing procedure of the texture analysis of 2nd embodiment. (A)-(D) is a figure which shows the image created by applying the texture analysis of 2nd embodiment, (E) and (F) is a figure which shows the MD image and FA image calculated by the conventional method. It is a figure which shows the simulation result of a diffusion model, (A) is a figure which shows FA of a some diffusion model, (B)-(E) is a figure which shows a texture index about the same model as (A).

<First embodiment>
First, an embodiment of an image processing apparatus will be described with reference to the drawings.
As shown in FIG. 1, the image processing apparatus 200 includes an image receiving unit 210 that receives image data acquired by an image capturing apparatus, a data analysis unit 220 that performs data analysis on the image data received by the image receiving unit 210, and a data analysis unit. The image generation unit 230 that generates a new calculation image that is the analysis result of the unit 220. The imaging device is, for example, an MRI device, and acquires a plurality of image data obtained by imaging under different imaging conditions, particularly spatial conditions such as an application axis of a gradient magnetic field. The image processing apparatus 200 further includes an input device 250 for inputting conditions and commands necessary for processing of each unit, and a display device (not shown) for displaying or storing data during or after processing. ) And a storage device 270 may be provided.

The data analysis unit 220 extracts the relationship (reciprocal position) of the spatial conditions for the plurality of image data with different spatial conditions received by the image reception unit 210, and based on the extracted mutual position, the voxel value of the image data. An index representing the distribution state of is calculated. Therefore, the data analysis unit 220 uses the relationship between the spatial conditions, that is, the mutual position calculation unit 221 that calculates the mutual position of each spatial condition, and the voxel value using the mutual position and each voxel value of the image data. And an index calculation unit 223 that calculates an index (texture index) that represents the distribution state of.

For example, if the space condition is represented by a vector such as a direction, the mutual position extracted by the mutual position calculation unit 221 includes the relative angle of each vector and the relative position of the coordinates of the base vector on the spherical surface. Is. The index calculation unit 223 selects, from the mutual positions obtained by the mutual position calculation unit 221, one or a plurality of vectors in which the mutual positions have a predetermined relationship, and uses each pixel value (physical quantity) in these vectors. Then, the feature amount of the texture analysis is calculated. Various texture feature amounts are known depending on the analysis method. For example, uniformity (ACM: Angular Second Moment), local difference degree (IDM: Inverse Difference Moment). The degree of local difference (Contrast), the degree of randomness (Entropy), the local correlation (Correlation), and the like are indexes (texture indexes) indicating the spatial distribution of physical quantities in each voxel. Therefore, the image generation unit 230 creates an image in which the index is the pixel value, and thus an image representing the characteristics of the spatial distribution is obtained.

The image generated by the image generation unit 230 may be displayed on the display device when the image processing device 200 includes a display device, or may be sent to the image pickup device (MRI device) 100 as image data, and is required here. According to the above, it may be displayed together with the original image data or a calculation image obtained by processing the original image data. The image data can be exchanged between the image processing apparatus 200 and the image capturing apparatus 100 via any known transfer means such as the Internet or a portable medium, and any of them may be adopted.

Further, the functions of the respective units constituting the image processing apparatus 200 can be realized by uploading software programmed in advance to a computer having a CPU, a GPU and a memory. Further, some functions may be realized by hardware such as ASIC and FPGA. When the image capturing apparatus 100 includes an image processing unit that creates an image from the acquired data, the image processing unit can also realize the above-described functions of the image processing apparatus 200.

FIG. 2A shows an outline of processing of the image processing apparatus of this embodiment. Further, as a comparative example, FIG. 2B shows an outline of the processing of the conventional method using the tissue model.

As shown in FIG. 2, regardless of the method of the present embodiment or the conventional method, first, a plurality of measurement values having different spatial conditions are input (S21, S31). The measurement value may be a signal measured by the image pickup apparatus 100 or image data obtained by reconstructing the signal, or image data obtained by further converting the signal into a predetermined physical quantity in the image pickup apparatus 100. In the former case, the image processing apparatus 200 processes the received image data to convert the signal into a physical quantity for each spatial condition, and creates a plurality of image data (S22, S32). Each voxel that constitutes each image data is a physical quantity in a specific spatial condition, for example, a physical quantity in a specific direction, and the purpose of the subsequent processing is to analyze the physical quantity in a plurality of spatial conditions and to calculate the spatial distribution of the physical quantity. Is to present the index or the image thereof (S25, S34). In the conventional method, a relational expression between a physical quantity and a direction is set using a model of spatial distribution assumed in advance (S33), and an index representing the spatial distribution of the physical quantity is calculated by approximating the relational expression. On the other hand, in this embodiment, first, the mutual spatial relationship (mutual position) is calculated only by the spatial condition (S23). Next, using the spatial relationship and the physical quantity in each direction, by applying the method of texture analysis, an index representing the characteristic of the spatial distribution of the physical quantity for each voxel, for example, uniformity, local difference degree, contrast, randomness. Degree, local correlation, etc. are calculated (S24). By imaging the index calculated for each voxel, it is possible to grasp whether the spatial distribution is uniform, random, or has a local correlation.

In the conventional method of presenting an index using a predetermined spatial distribution model as shown in FIG. 2B, the spatial distribution model assumed by the relational expression between the physical quantity and the directionality is limited. Therefore, if the target tissue structure and the spatial distribution of the physical quantity of the assumed model do not match, the diagnostic interpretation of the presented index will be difficult. For example, in the case of assuming a model in which a physical quantity follows an elliptic three-dimensional distribution, if the physical quantity actually has a spatial distribution that is a combination of two elliptical three-dimensional distributions, then the hypothetical model uses a space for the physical quantity. The characteristics of the distribution cannot be expressed accurately. On the other hand, according to the image processing apparatus of the present embodiment, since the characteristics of the spatial distribution of physical quantities are obtained without using a model or approximation, information useful for diagnosis is presented regardless of the target tissue. be able to.

<Second embodiment>
An embodiment in which the image pickup apparatus 100 is an MRI apparatus and the image processing unit of the MRI apparatus has the function of the image processing unit 200 will be described based on the above-described first embodiment.

First, the overall configuration of the MRI apparatus to which this embodiment is applied will be described.
As shown in FIG. 3, the MRI apparatus 100 includes, as a measuring unit 10, a static magnetic field magnet 101 that generates a static magnetic field, a gradient magnetic field coil 102 that generates a gradient magnetic field that gives a magnetic field gradient to the static magnetic field, and a static magnetic field space. Generated from the subject 103, and an RF transmission coil 107 for applying a high-frequency magnetic field pulse that causes nuclear magnetic resonance to nuclear spins (hereinafter, simply referred to as spins) of atoms contained in the tissue constituting the subject 103 arranged in An RF receiving coil 108 for receiving a nuclear magnetic resonance signal and a table 115 for loading the subject 103 into the static magnetic field space are provided. The measuring unit 10 may include a shim coil 113 for improving the uniformity of the static magnetic field.

Although the gradient magnetic field coil 102 is omitted and shown as one block in the figure, it is configured by coils in three axial directions orthogonal to each other, and the gradient magnetic field coils of the respective axes are connected to the individual gradient magnetic field power supplies 105. By controlling the drive current supplied from each gradient magnetic field power source, a gradient magnetic field of a desired direction (application direction) and a desired magnitude can be generated. The shim coil 113 is connected to the shim power supply 114 and generates a predetermined gradient magnetic field by the shim current.

The RF transmission coil 107 is connected to the transmission unit 106 including an oscillator that generates an RF signal of a magnetic resonance frequency, a high frequency amplifier, and the like, and generates an RF pulse having a predetermined shape. The RF receiving coil 108 is connected to a receiving unit 109 including a high frequency amplifier, a quadrature detection circuit, an A/D converter, etc., and sends the received nuclear magnetic resonance signal to the receiving unit 109.

The MRI apparatus 100 includes a sequencer 104 that controls the operations of the gradient magnetic field power source 105, the transmission unit 106, the reception unit 106, and the shim power source 114 described above according to a predetermined pulse sequence, and a computer 110 that controls the entire apparatus and performs signal processing. Prepare In the present embodiment, a pulse sequence for diffusion weighted imaging described below is executed as the pulse sequence.

The computer 110 includes a processor such as a CPU and a GPU and a memory, and further, as an external device, an input device 111 for a user to input commands and conditions necessary for device control, signal processing, and the like, processing results, and the like. Is provided with a display device 116 for displaying, a storage device 112 for storing data necessary for control, a program, data being processed, and the like as necessary.

The computer 110 uploads pre-programmed software to a processor and executes the software to perform calculations for various controls and image processing. The function executed by the computer 110 specifically includes a control unit 30 (measurement control unit 31, display control unit 33) and an image processing unit 20, as shown in the block diagram of FIG. Includes an image reconstruction unit 21, an image analysis unit 22, and an image generation unit 23.

The measurement control unit 31 sends the imaging sequence (pulse sequence) and imaging conditions set by the user to the sequencer 104 via the input device 111, and controls the operation of the measurement unit.

The image reconstructing unit 21 reconstructs an image by performing calculations such as Fourier transform on the nuclear magnetic resonance signals collected by the receiving unit 106. The image analysis unit 22 has a function of performing analysis on a reconstructed image and creating a new image, and a physical quantity calculation unit 227 for calculating a physical quantity of a subject tissue based on the reconstructed image, and FIG. The image processing apparatus 200 includes the same functional unit as the image analysis unit 220, that is, the mutual position calculation unit 221 and the index calculation unit 223. In the case of diffusion weighted imaging, the physical quantity calculated by the physical quantity calculation unit 227 includes various quantities such as an apparent diffusion coefficient (ADC) and a diffusion kurtosis coefficient (K). Specific processing of these image processing units 20 will be described later.

The image generation unit 23 uses the index calculated by the index calculation unit 223 to generate a new image having the index as a voxel value.
The display control unit 33 controls to display the images generated by the image reconstruction unit 21, the physical quantity calculation unit 227, and the image generation unit 23 on the display device 116 in a predetermined display form together with information about the subject.

Next, the operation of the MRI apparatus 100 will be described based on the above configuration. Here, the operation of the image processing unit 20 will be mainly described by taking the case of performing diffusion-weighted imaging as an example.
As shown in FIG. 5, the operation flow includes an imaging step S51, an image reconstruction step S52, a physical quantity calculation step S53, a mutual position calculation step S54, a texture index calculation step S55, and an image generation step S56. The details of each step will be described below.

[Imaging step S51]
Under the control of the sequencer 104, the measurement unit 10 executes a pulse sequence used for diffusion weighted imaging, and acquires a diffusion weighted image. An example of the diffusion weighted pulse sequence is shown in FIG. The illustrated pulse sequence is a spin echo type EPI sequence, in which an RF pulse 602 for exciting a spin is first applied together with a slice selection gradient magnetic field 601 and then an RF pulse 605 for reversing the spin is applied together with a slice selection gradient magnetic field 604. Between the excitation RF pulse 602 and the inversion RF pulse and after the application of the inversion RF pulse, gradient magnetic field pulses (MPG pulses) 603 and 606 with high intensity are applied. Then, after applying the phase-encoding gradient magnetic field and the read-out gradient magnetic field dephasing pulses 607 and 608, respectively, while inverting the intensity of the read-out gradient magnetic field pulse 612 and applying the phase-encode gradient magnetic field pulse 611 in a blip shape, The nuclear magnetic resonance signal is collected as the echo signal 610. One excitation may collect all the echo signals necessary for one image, or may be divided into a plurality of excitations and collected. The echo signal is arranged at the position of the k-space data determined by the phase encoding gradient magnetic field and the readout gradient magnetic field at the time of acquisition, and one piece of k-space data corresponding to one image is collected.

In FIG. 6, MPG pulses are applied simultaneously on three axes (Gz, Gy, Gx), but different axes are applied to collect echo signals a plurality of times, and a plurality of k-space data with different MPG pulse conditions are acquired. collect. The number of MPG pulse application axes is not particularly limited, but is usually about 6 to 30, and imaging may be performed with one application axis having different b values of 2 or more. The b value (s/mm ² ) is selected from values such as 500, 1000, 1500, 2000, 2500. The pulse sequence in which no MPG pulse is applied (b value=0) is common to each application axis and may be executed once. The reconstructed image of k-space data obtained by imaging this pulse sequence is referred to for each application axis. Let it be an image.

[Image reconstruction step S52 and physical quantity calculation step S53]
The image reconstruction unit 21 performs an operation such as Fourier transform on the plurality of k-space data collected in step S51 to create image data (S52). That is, the image data is composed of a plurality of image data that differ depending on the application axis of the MPG pulse and the set b value.

式（１）中、ｂはｂ値、Ｓ_０はＭＰＧパルスを印加しない参照画像の信号値、Ｄ_ｎ,appはみかけの拡散係数である。 The physical quantity calculation unit 227 calculates a physical quantity such as an apparent diffusion coefficient using a plurality of image data (S53).
First, assuming that there are N application axes, the signal value S(n,b) of the nth (any one of 1 to N) application axis is given by the following equation (1).

In Expression (1), b is the b value, S ₀ is the signal value of the reference image to which the MPG pulse is not applied, and D _n,app is the apparent diffusion coefficient.

Further, the physical quantity calculation unit 227 may use the diffusion coefficient for each axis to calculate the average diffusion coefficient MD and the diffusion anisotropy FA as in the conventional method. These calculations are not mandatory.
If a three-dimensional spheroidal distribution model is used here, D _n can be expressed by the equation (2) using the diffusion tensor D _i,j , and then the average diffusion coefficient MD can be calculated using the equation (3). The diffusion anisotropy FA can be calculated by using the equation (4).

式（３）中、λは拡散テンソルの固有値であり、添え字ｉは３軸方向のいずれかを意味する（以下、同じ）。

In Expression (3), λ is the eigenvalue of the diffusion tensor, and the subscript i means any of the three axial directions (hereinafter the same).

The image generation unit 23 generates image data having the physical quantity calculated by the physical quantity calculation unit 227 as a pixel value. The physical quantity itself is information representing, for example, a diffusion distribution state of cerebral blood flow, and serves as a diagnostic index for acute cerebral infarction or white matter lesion. For example, since the white matter nerve fibers have a tubular structure, the diffusion increases along the traveling direction, so that the white matter can be emphasized in a highly anisotropic region. However, as the anisotropy becomes smaller, the distribution shape tends to become more complicated, and the index of anisotropy cannot sufficiently index the diffusion distribution. In the next steps S54 to S55, the diffusion distribution, which cannot be indexed by these physical quantities, is indexed by the feature quantity of the texture analysis.

[Mutual Position Calculation Step S54]
When the image analysis unit 22 receives the physical quantity image data (for example, the diffusion coefficient image) generated in step S53 described above, as a first step of image analysis, the mutual position calculation unit 221 determines that the relative positions of the plurality of application axes with respect to each other. Calculate (relationship) and convert to data. The information about the application axis of the MPG pulse received by the image processing unit 23 is, for example, (1, 0, 0), (0, 0) as the intensity (Gx, Gy, Gz) of each axis Gx, Gy, Gz of the gradient magnetic field pulse. 1,0), (0,0,1) or (0.525,0,-0.851), (-0.525,0,-0.851), (0.851,0.525,0) ), etc., and their relationship is unknown.

The mutual position calculation unit 221 calculates relative positions such as distances and angles with each other using the information of the applied axes. Several methods can be used to calculate the relative position. Here, a method using Delaunay triangulation as an example of the relative position will be described with reference to FIG. 7.

First, the MPG application axis direction (diffusion direction (ex, ey, ez)) is standardized to form a vector of magnitude 1. Note that the application axis direction of the MPG pulse can be rephrased as the diffusion direction if it is defined in association with the physical quantity called the diffusion coefficient, and is therefore referred to as the diffusion direction below. As a result, the MPG pulse application axis direction, that is, the diffusion direction can be plotted on a sphere having a radius of 1 as shown in FIG. Next, a triangular surface mesh is created from the normalized coordinates in the diffusion direction based on the Delaunay triangulation conditions (FIG. 7(B)). The Delaunay condition is that the circumscribed circle of each triangle does not include other points. In this triangular surface mesh, a plurality of adjacent diffusion directions are connected to one diffusion direction as one side of the triangle (FIG. 7C).

[Texture index calculation step S55]
Next, we quantify the diffusion information by applying the texture analysis. The texture analysis method uses a gray level co-occurrence matrix, a method of calculating a feature amount using a density histogram, and a histogram feature amount using a difference in density value between a pair of pixels having a predetermined relative positional relationship. There is a method of calculating a feature amount (a method of using a difference statistic) similarly to the above, but an example in which a gray level co-occurrence matrix (GLCM: gray level cooccurrence matrix) is applied will be described. First, as a reference, a texture analysis using a conventional GLCM will be briefly described with reference to FIG.

It is assumed that there is image data 800 as shown in FIG. Small squares in the quadrangle shown on the left side of FIG. 8 correspond to pixels (voxels), the numbers therein are pixel values, and the pixel values are standardized to 1 to 5. The GLCM 801 shown on the right side of FIG. 8 is a matrix (matrix size: 5×5) in which the vertical and horizontal directions correspond to the standardized pixel values, the vertical direction is the left-to-right relationship of the pixels of the image data 800, and the horizontal direction is from the right side of the image. The relationship on the left is shown. Of the pixels of the image data 800, the pixel having a pixel value of 1 is focused on, and the combination in which the value of the pixel adjacent to the right thereof is 5 is counted. In the figure, a combination surrounded by a dotted line corresponds to the combination. The count number is set to a value of 1 in the row and 5 in the column of the GLCM 801. In the same manner, a combination (pixel pair) of all pixels having a row and column relationship of the GLCM is found in the same manner, and the count number is set to the value of the corresponding row and column of the GLCM to complete the GLCM. Using the matrix thus obtained, the texture feature amount is calculated.

The texture analysis adopted in this embodiment is an application of this GLCM method, but the general texture analysis applied to medical images and the like is to obtain the feature amount of a two-dimensional image. On the other hand, here, the feature amount of the spatial distribution of the physical amount is calculated. That is, the target data is not the image data itself obtained by imaging, but a physical quantity image representing a physical quantity such as a diffusion distribution obtained by calculation from the image data itself. Further, the positional relationship of the pixels referred to by the GLCM is not the relative position of the pixel in the two-dimensional coordinate of the image data as in the conventional texture analysis, but the diffusion direction connected by the triangular surface mesh calculated by the mutual position calculation unit 221. It is a positional relationship.

Hereinafter, the specific procedure of the texture analysis of the present embodiment will be described by taking the case where the physical quantity image is a diffusion coefficient image as an example. FIG. 9 shows a flow showing the procedure, and FIG. 10 shows an explanatory view.

First, the diffusion coefficient in each application axis direction is standardized (S91). Specifically, the maximum value of the measured (calculated) diffusion coefficient is set to 1. Next, the matrix size of GLCM is determined (S92). The matrix size corresponds to the number of gradations (for example, M) when the diffusion coefficient is expressed in gradations, and is, for example, about 10 to 50. The matrix size may be set in advance by a program that executes the mutual position calculation algorithm, or may be selectable by the user. Further, since the step size is determined by the matrix size of the GLCM, the diffusion coefficient can be calculated as a fixed value without being standardized.

Next, as shown in FIG. 10, in the triangular surface mesh, one or a plurality of diffusion directions in the vicinity of one diffusion direction are selected (S93), and GLCMs are created respectively (S94). In the example shown in FIG. 10, with respect to the diffusion direction p0, a plurality of (here, six) diffusion directions n1 to n6 are selected from the closest diffusion direction to the sixth closest diffusion direction. In the triangular surface mesh, triangles (six triangles in the figure) having the vertex in the diffusion direction p0 are determined, and therefore the diffusion directions at the vertex positions other than the diffusion direction p0 of these triangles may be selected. .. When the number of neighboring diffusion directions to be selected is predetermined, a predetermined number of neighboring diffusion directions may be selected based on the distance (distance from the diffusion direction p0 or mutual distance) from the selected neighboring diffusion directions. ..

Details of GLCM creation (S94) for each voxel will be described with reference to the flow of FIG. The process is sequentially executed for each of a plurality of voxels (for example, 256 voxels, 512 voxels) included in the image. In the processing of one voxel, the processing is sequentially executed for each of the plurality of application axis directions.

First, if a plurality of (here, six) neighboring diffusion directions are selected for one diffusion direction (S941), GLCM1 to GLCM6 representing the relationship with the pixel value of the pixel in the diffusion direction p0 are created (S942). .. The vertical axis and the horizontal axis of these GLCMs respectively correspond to the gradation of the diffusion coefficient, as described above. Then, for example, in the case of GLCM1, the pixel values (standardized value of the diffusion coefficient: standardized pixel value) of the diffusion coefficient images in the diffusion direction p0 and the diffusion direction p1 are compared between the same voxels (k-th voxel). , Which position of the GLCM1 the combination (pixel value pair) is, and counts the determined GLCM1 position. With respect to the same voxel, the same process is performed for the diffusion direction p0 and the other diffusion directions p2 to p6, and the Cowout numbers are input to the GLCMs 2 to 6.

The above processes S941 and S942 are repeated by changing the diffusion direction p0 to another diffusion direction, six neighboring diffusion directions are selected in the same manner as the diffusion direction p0, and pixel value pairs are set for the same voxel (kth voxel). Determines which position of GLCM is 1 to 6 and puts the count value. Finally, similar processing is performed for all diffusion directions p0 to pN (N is the number of application axes), and GLCMs 1 to 6 having count values as elements are completed.
Next, these GLCMs are averaged (S943), and an average GLCM is created as shown on the left side of FIG. This completes the creation of the average GLCM for one voxel.

The processes S941 to S943 described above are executed by changing voxels, and finally the average GLCM is created for all voxels (S944). Note that the process for each voxel may be performed for all voxels forming the diffusion coefficient image in each direction, but, for example, when the target region is a specific part such as the brain, threshold processing is performed, and the signal value is The processing of the pixels below the threshold value or the background pixels may be omitted. Alternatively, a mask for selecting only the target area from the diffusion coefficient image may be applied, and only the pixels in the selected area may be processed.

Finally, the feature amount is calculated for each voxel by using the average GLCM matrix p(i,j) (FIG. 9: S95). As the feature amount, a known feature amount for texture analysis can be obtained. Specifically, there are uniformity (ASM), local difference degree (IDM), contrast (Contrast), randomness (Entropy), local correlation (Correlation), etc., and the following equations (5) to (9) are used. It can be calculated. In these equations, i and j are the coordinates of each pixel forming the pixel pair (i,j), M is the GLCM matrix size, and p(i,j) is the average GLCM matrix.

式（９）において、μはi, jの平均値、σはi,jの標準偏差である。

In Expression (9), μ is the average value of i and j, and σ is the standard deviation of i and j.

Through the above processing S94 and S95, the analysis by the image analysis unit 22 and the calculation of the index (FIG. 5: S55) are completed.

[Image generation step S56]
The image generation unit 23 creates a texture image having the pixel value of the characteristic amount (index) calculated by the image analysis unit 22. The display controller 33 converts the texture image generated by the image generator 23 into display image data, and displays it on the display device 116 together with other supplementary information and the like. At this time, when the diffusion coefficient image and the FA image are created in step S53, these images can be displayed in parallel or selectively.

The texture image displayed in this way is an image that shows the spatial characteristics of the diffusion distribution, such as the uniformity and the degree of randomness of the diffusion distribution. Also, it is possible to provide information that is easy to understand.

[Effects of First Embodiment]
Quantitative images of ASM, IDM, Entropy, and Correlation created using the texture analysis of the present embodiment are shown in FIGS. 12(A) to (D). These images were created from diffusion-weighted images obtained by imaging with 21 MPG axes and 1 b value (1000 [s/mm ² ]). For reference, MD (average diffusion coefficient) images and FA images created from the same diffusion-weighted image are shown in FIGS. 12(E) and 12(F). In the MD image, since the diffusion coefficients are averaged, only the intensity of diffusion can be known, and direction information cannot be obtained. Further, in the FA image, the white matter with a clear diffusion direction is clearly drawn, but the contrast is reduced at the portion where the plurality of strong diffusion directions shown by the arrows intersect (the area of the intersecting fibers), and clear information cannot be obtained. On the other hand, in the images of FIGS. 12A to 12D, for example, by comparing the ASM image and the IDM image or by observing the image of the degree of disorder, whether the diffusion distribution is uniform or scattered. It is possible to obtain information such as, and from the local correlation image, information such as whether there is a correlation between the diffusion of a certain region and the diffusion of an adjacent region.

FIG. 13 is a diagram showing a result of confirming the superiority of the quantitative value image calculated in this embodiment by simulation. In the simulation, various diffusion models with different diffusion distributions were created, and among them, a plurality of models in which the average diffusion coefficient and the diffusion anisotropy were almost uniform were selected. FIG. 13A is a graph showing the diffusion anisotropy of the selected model, and the horizontal axis is the model number. With respect to the plurality of diffusion models, 21-axis diffusion directions were set, and texture analysis using relative positions in the diffusion directions was performed by the method of the first embodiment to calculate ASM, IDM, randomness, and local correlation. 13B to 13E are graphs showing the results. Also in these graphs, the horizontal axis is the model number. As can be seen from these results, even if the diffusion models of FA have almost the same values, the difference in diffusion distribution appears as a clear difference by adopting the texture analysis method of the present embodiment. I understand.

The embodiment in which the image processing of the present invention is realized by the image processing unit of the MRI apparatus has been described above. However, in the MRI apparatus, only the diffusion-weighted image is captured, or the capturing to the physical quantity calculation process is performed, and the subsequent processes are performed. It may be executed by an image processing apparatus different from the MRI apparatus. In that case, for example, the function of the image generation unit 23 and the display on the display device may be executed by the MRI apparatus in parallel with the image processing apparatus.

According to the present embodiment, the diffusion image is analyzed by using the spatial relationship of the MPG axis, that is, the diffusion direction, instead of using the conventional model assuming a predetermined tissue, and thus the diffusion by the tissue is performed. Even if the coefficient distribution is not uniform and a complicated portion or organization is used, the complexity can be shown as an index, and detailed information that cannot be obtained by the conventional FA or MD can be obtained. In addition, it is possible to present information for grasping the distribution state of diffusion calculation without being restricted by models and approximations.

Further, the triangular surface mesh adopted in the present embodiment can uniquely determine a neighboring diffusion direction having a predetermined spatial relationship with one diffusion direction (for example, in the triangular surface mesh, the diffusion direction becomes the apex). Since a plurality of triangles are defined, the diffusion directions of the other vertices of these triangles are uniquely determined), and it becomes possible to easily select the diffusion direction in the subsequent texture analysis processing.

<Modification of physical quantity data>
In the first embodiment, the physical quantity data to be processed by the image processing unit 20 is the diffusion coefficient (ADC) of each application axis, but a kurtosis coefficient or other data may be used as the physical quantity data.

For example, the equation (1) representing the attenuation of the signal by the b value is a mathematical model assuming that the diffusion has a Gaussian distribution, but can be represented by the equation (10) assuming a non-Gaussian distribution.

式中、Ｋnは拡散分布の尖度を表す尖度係数である。１つのＭＰＧ印加軸について、ｂ値が異なる複数の条件で信号を取得することにより、式（１０）から、拡散係数Ｄnと尖度係数Knを算出することができる。こうして算出される尖度係数も拡散係数の場合と同様に、その空間分布をテクスチャー解析によって提示することができる。

In the equation, Kn is a kurtosis coefficient representing the kurtosis of the diffusion distribution. The diffusion coefficient Dn and the kurtosis coefficient Kn can be calculated from Expression (10) by acquiring the signals under a plurality of conditions with different b values for one MPG application axis. The kurtosis coefficient thus calculated can also be presented by a texture analysis, as in the case of the diffusion coefficient.

In addition, there is a mathematical model showing the attenuation of a signal in which two diffusion coefficients are considered (for example, Expression (11)). The fast diffusion coefficient (D ^* ) and the slow diffusion coefficient (D) defined by these equations can be processed as physical quantity data in the same manner as the diffusion coefficient, and the index of the spatial distribution thereof can be calculated.

<Modified example of mutual position calculation>
In the first embodiment, the case where the mutual position calculation unit 221 calculates the triangular surface mesh as the data representing the spatial relationship of the application axis of the MPG pulse has been described, but the spatial relationship is not limited to the triangular surface mesh. , As long as it is a data representation of each relationship, it may be information on the angle between the application axes.

In this case, the steps of mutual position calculation and index calculation (S54, S55) are changed as follows. First, in the mutual position calculation step S54, the MPG application direction is standardized to the size 1 as in the first embodiment. Next, the angle formed by each application direction is calculated. The application direction is determined by the combination (ax, by, cz) of the application strengths of the gradient magnetic fields Gx, Gy, Gz. For example, if the application direction is the Gx direction, a=1 and b=c=0 (1,0 , 0), if the application direction is the Gy direction, (0, 1, 0) is expressed as a vector, and if it has an angle with respect to the three orthogonal axes, (0.525, 0, −0.851). , (-0.525,0,-0.851), (0.851,0.525,0), etc. The angle can be obtained by normalizing the magnitudes of these vectors and then taking the inner product.

Next, in the index calculation step S55, when selecting the neighborhood application directions of the respective application directions (diffusion directions), a plurality of neighborhood application directions are selected in ascending order of angle. Once the application direction is determined, GLCM is created for each application direction, averaged, and the quantitative value of the texture analysis is calculated from the average GLCM as in the first embodiment.

<Modified example of texture analysis>
In the first embodiment, the method of using the gray level co-occurrence matrix (GLCM) has been described as the texture analysis method. However, the length of consecutive same pixel values is counted, and GLSZM(Gray) having the number as a matrix element is used. A matrix such as Level Size Zone Matrix) may be used.

Further, as a method of texture analysis, a feature amount is calculated using a density histogram, or a feature amount is calculated in the same manner as the histogram feature amount by using a difference in density value between pixel pairs having a predetermined relative positional relationship. There is a method (a method of using a difference statistic) or the like, which may be applied. As an example, a modification to which the method using the difference statistic is applied will be described.

Also in this modification, the mutual position calculation unit 221 calculates the mutual position of the MPG application axes represented by a triangular surface mesh or the like. Next, as in the first embodiment, the index calculation unit 223 selects a plurality of neighboring diffusion directions for one diffusion direction, and for each of them, a difference (density) between pixel values of a corresponding pixel pair (pixel pair). The difference) is calculated, and the pseudo probability P(l) of the difference is calculated. The probability P(l) is used to calculate a feature amount such as an average value (MEN), a variance (VAR), a contrast (CNT), a skewness (SKW), and a kurtosis (KRT) from a density histogram. The feature amount based on the difference statistic is calculated by using the equation (1).

Although the second embodiment and the modified examples thereof have been described above, these modified examples can be appropriately combined, and such combinations are also included in the modified examples of the present embodiment.

<Third embodiment>
In the second embodiment, the case where the image acquired by the MRI apparatus is calculated from the diffusion weighted image is explained. However, the image targeted by the image processing apparatus of the present invention or the image processing unit of the MRI apparatus is a diffusion weighted image. Not limited.

For example, a T2* image having a lateral relaxation time as a pixel value and a magnetic susceptibility image are both considered to have a spatial distribution due to the dependence on the direction of the static magnetic field. Therefore, a plurality of images captured by changing the angle of the subject (for example, the angle of the head) with respect to the static magnetic field are images of physical quantities with different spatial conditions.

Therefore, similar to the first and second embodiments, the relative positions of a plurality of positions (here, the angle with respect to the static magnetic field) are converted into data, and the pixel pairs between the images at the adjacent positions are analyzed by the texture analysis method. By doing so, it is possible to image the features of the spatial distribution of T2* and magnetic susceptibility.

Although the embodiment of the present invention and the modified example thereof have been described above, one feature of the present invention is that the measured physical quantity and the spatial information inherent therein are modeled and analyzed on the premise of their relational expressions. Instead, first calculate the relationship (relative position) of the spatial information itself and use that information to analyze the predetermined spatial information and the associated physical quantity by applying a texture analysis method. However, the processing is not limited to the above embodiment as long as the processing has this feature. Further, the feature of the present invention is not to perform texture analysis using the position of the pixel pair in the two-dimensional image data, but to use a three-dimensional relative position and create a two-dimensional histogram for the pixel pair having a predetermined relative position. There is also the point of using a new analysis method. As a result, it is possible to present an index indicating the spatial distribution state of a physical quantity represented by a tensor that cannot be described in two-dimensional coordinates.

100: imaging device (MRI device), 10: measurement unit, 20: image processing unit, 21: image reconstruction unit, 22: image analysis unit, 23: image generation unit, 30: control unit, 31: measurement control unit, 33: display control unit, 110: computer, 111: input device, 112: storage device, 116: display device, 200: image processing device, 220: image analysis unit, 221: mutual position calculation unit, 223: spatial index calculation unit 227: Quantitative value calculation unit

Claims (14)

A receiving unit that receives a plurality of image data in which the physical quantity is a pixel value and the physical quantity has spatial information,
A mutual position calculation unit that extracts the relationship of spatial information in each of the plurality of image data,
An image processing apparatus, comprising: a texture analysis unit that analyzes the distribution of the physical quantity and calculates a texture index using the spatial relationship extracted by the mutual position calculation unit. The image processing apparatus according to claim 1, wherein
The image processing apparatus is characterized in that the plurality of image data are a plurality of image data created from diffusion-weighted images captured by the magnetic resonance imaging apparatus with different application axes of MPG pulses. The image processing apparatus according to claim 2, wherein
The image processing device, wherein the spatial information is the direction of the application axis of an MPG pulse used for diffusion weighted imaging. The image processing apparatus according to claim 2, wherein
The physical quantity is any one of coefficients related to diffusion in a one-dimensional diffusion signal strength model including an apparent diffusion coefficient (ADC), a kurtosis coefficient (K), a slow diffusion coefficient, and a fast diffusion coefficient. Image processing device. The image processing apparatus according to claim 1, wherein
The image processing apparatus, wherein the plurality of image data are images captured by a magnetic resonance imaging apparatus at a subject position having a plurality of different angles with respect to a static magnetic field direction. The image processing apparatus according to claim 1, wherein
The mutual position calculation unit performs Delaunay triangulation using the spatial information, and the texture analysis unit uses position information connected by Delaunay triangulation as the spatial relationship. apparatus. The image processing apparatus according to claim 2, wherein
The image processing apparatus, wherein the mutual position calculation unit uses a vector representing an application axis direction of the MPG pulse as the spatial information, and calculates an inner product of the vectors as the spatial relationship. The image processing apparatus according to claim 1, wherein
The texture analysis unit selects, for each voxel forming the image data, one or a plurality of neighboring spatial information in the vicinity of each spatial information by using the spatial relationship, and selects the plurality of selected spatial information. An image processing apparatus, characterized in that each physical quantity is analyzed using a gray level co-occurrence matrix to calculate the texture index. The image processing apparatus according to claim 8, wherein
The image processing apparatus, wherein the texture analysis unit normalizes the physical quantity according to a matrix size of the gray level co-occurrence matrix. The image processing apparatus according to claim 8, wherein
The image processing apparatus, wherein the texture index includes any one of uniformity, local difference, contrast, randomness, and local correlation. The image processing apparatus according to claim 1, wherein
The image processing device further comprising an image generation unit that generates an image having a pixel value of a texture index calculated for each voxel by the texture analysis unit. The image processing apparatus according to claim 1, wherein
The reception unit receives, as the plurality of image data, a plurality of diffusion-weighted images captured by the magnetic resonance imaging apparatus with different application axes of MPG pulses,
The image processing apparatus further comprising a physical quantity calculation unit that creates image data having the physical quantity as a pixel value from the diffusion weighted image. A reception unit that receives a diffusion-weighted image captured using a plurality of MPG pulses with different application axes using the magnetic resonance imaging apparatus or a plurality of image data created from the diffusion-weighted image,
A mutual position calculation unit that calculates a mutual relationship in the application axis direction of the plurality of MPG pulses as a spatial relationship,
For each voxel of the image, one or more application axis directions having a predetermined relationship with each of the plurality of application axis directions are selected, and the texture index is calculated using the pixel value for the selected application axis direction. An image processing apparatus, comprising: a texture analysis unit that uses the texture index as a voxel pixel value. A measuring unit that measures a nuclear magnetic resonance signal using a plurality of MPGs having different application axes;
An image reconstruction unit that creates an image using the nuclear magnetic resonance signal measured by the measurement unit,
An image processing unit for processing the image created by the image reconstruction unit;
Equipped with
The magnetic resonance imaging apparatus, wherein the image processing unit includes the image processing apparatus according to claim 1.

Images