JP2005253702A - Image reconstruction method, apparatus and program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To remove image components other than a focal plane superimposed on the focal plane image. <P>SOLUTION: A processing part inputs an NMR signal or the like imaged by superimposing information of respective cross sections of a subject by a nuclear magnetic resonance imaging method, from a storage part or the like (S100). The processing part finds a first focusing image using a desired depth coordinate as a focal plane and its absolute value amplitude image based on the NMR signal (S301 and 102). The processing part finds a point image distribution function indicating an effect which the image components other than the focal plane superimposes on the focusing image and its absolute value image (S302). The processing part executes reverse filtering based on the absolute value image of the point image distribution function and the absolute value amplitude image of the focusing image to find an incoherent image formation approximate image (hereinafter referred to as approximate image) with improved image quality (S303). The processing part finds a modified function based on a ratio of the second focusing image to the first focusing image which are obtained based on the approximate image and the point image distribution function, and corrects the approximate image according to the modified function to find a reconstruction image removed with the images other than the focal plane (S305). <P>COPYRIGHT: (C)2005,JPO&NCIPI

Description

  The present invention relates to an image reconstruction method and apparatus, and an image reconstruction program, and in particular, one or several images reconstructed from signals in magnetic resonance imaging using Fresnel diffraction and other image sensing methods. The present invention relates to an image reconstruction method and apparatus in magnetic resonance imaging using a three-dimensional diffraction distribution, and an image reconstruction program, in which image components other than the focal plane are removed from (focused image) and only the focal plane image component is extracted.

  In order to acquire a three-dimensional image in magnetic resonance imaging (MRI, Magnetic Resonance Imaging), usually three-dimensional data collection is performed. Therefore, generally a long imaging time is required. On the other hand, in a magnetic resonance imaging method (hereinafter referred to as Fresnel diffraction type imaging method) that obtains a signal having the same shape as the three-dimensional diffraction distribution or uses the Fresnel diffraction method, even when the image is reconstructed from a two-dimensionally scanned signal. An image focused on an arbitrary depth depth can be obtained according to the perspective parameter to be used. This is similar to the principle that an image focused at an arbitrary position can be obtained by adjusting the focal length of a lens with a camera or the like. If three-dimensional information can be obtained by two-dimensional data acquisition scanning, it is possible to reduce MRI imaging time, which takes a long time, to reduce physiological effects by RF (Radio Frequency), and to reduce the mental burden on the subject. Its utility is extremely large.

  Further, in an apparatus for obtaining an image of a surface of a subject at an arbitrary depth by X-rays, the subject is irradiated with X-rays and detected by a light receiving element functioning as a means for imaging X-rays transmitted through the subject. An X-ray imaging apparatus that obtains an image of a surface at an arbitrary depth is disclosed (for example, see Patent Document 1). Japanese Patent Application Laid-Open No. H10-228707 describes that a distance adjustment unit adjusts the distance between a light collecting unit and a subject to obtain an image of a surface at an arbitrary depth of the subject. Further, in Patent Document 1, in order to extract a stereoscopic image, data is extracted while relatively moving a tomographic plane parallel to the xy plane along the z axis at a certain interval, and the plane data is continuously or It is described that a three-dimensional image can be directly measured without depending on calculation by intermittently measuring and superimposing on a computer.

Further, the literature of the present inventors describes an NMR (Nuclear Magnetic Resonance) Fresnel diffraction imaging method capable of acquiring an image having information in the depth direction by collecting two-dimensional data (for example, Non-Patent Document 1). reference). In this method, by devising the shape of the magnetic field applied at the time of MRI measurement, a signal having the same shape as the three-dimensional diffraction distribution of light can be obtained, and the holography principle is applied to focus on an arbitrary depth from a signal scanned two-dimensionally. A combined image can be obtained, but it is described that a blurred image out of focus is superimposed on the obtained image.
Further, there is a document describing an iterative formula for deconvolution integral based on Bayesian theory (see Non-Patent Document 2).
JP 2000-135211 A Satoshi Ito, Yukazu Ono, Yoshiaki Uemura, Yoshifumi Yamada, “Image Reconstruction by NMR Fresnel Diffraction Imaging”, MEDICAL IMAGEING TECHNOLOGY vol. 18, no. 6, November 2000, pp. 817-827 F. Tsumuraya, “Deconvolution based on the Wiener-Lucy chain algorithm: an approach to recover local information losses in the deconvolution procedure”, J.Opt.Soc.Am.A, Vol.13, No.7, pp.1532- 1536 (1996).

  In order to acquire a three-dimensional image, usually three-dimensional data collection is performed, so that the imaging time is generally increased. For example, as in the apparatus described in Patent Document 1, in order to obtain a three-dimensional image of a subject from a plurality of two-dimensional images, it takes a measurement time to acquire each two-dimensional image.

  Because NMR imaging is so-called “perspective” that captures information inside an object, light is emitted not only from the surface but also from the inside of the object, as is the case with natural objects handled by optical holography. The format is Therefore, the information of each cross section in the depth direction is overlapped and imaged. That is, although the focal plane image is clear in the reconstructed image, there is a problem that the cross-sectional image components before and after the superimposition are superimposed in a blurred shape according to the distance from the focal plane position. For example, in image diagnosis, the presence of such a background image may make it difficult to distinguish between a normal part and an abnormal part. For this reason, it is desirable to extract only the focal plane image as much as possible and to remove image information components other than the focal plane.

  In view of the above, the present invention provides an image reconstruction method and an image reconstruction program for removing a background image component superimposed on a focal plane image from a focal plane image group obtained in the Fresnel diffraction type imaging method. Objective. It is another object of the present invention to obtain a clearer image by removing image information components other than the focal plane, and to easily identify, for example, a normal site and an abnormal site in image diagnosis.

  Another object of the present invention is to significantly reduce the MRI measurement time. More specifically, it is possible to perform three-dimensional information by two-dimensional data collection scanning, and to reduce MRI imaging time that takes a long time. As a result, the present invention also aims to improve throughput due to the increase in the number of MRI examinees per day, reduce the cost of consultation, and alleviate mental distress of subjects placed in an obstructive environment. is there. Another object of the present invention is to reduce the physiological effects of RF (neural stimulation due to eddy current generation in the living body, increase in living body temperature) on the subject.

  The present invention is an application of the principle of reproducing a stereoscopic image from a hologram in which a wavefront is two-dimensionally recorded in optical holography to an NMR imaging method. The present invention removes the background image component superimposed on the focal plane image from the focal plane image group obtained by the Fresnel diffraction type imaging method by inverse filtering of incoherent imaging approximation and the maximum likelihood method.

  Specifically, first, a focal plane image with improved image quality is obtained by performing inverse filtering of incoherent imaging approximation on a focal plane image on which blur obtained by the Fresnel diffraction type imaging method is superimposed. Next, the focal plane image without blur is obtained by applying the maximum likelihood method using this image as an initial value. The convergence time of the maximum likelihood method can be shortened by using an image obtained by inverse filtering of incoherent imaging approximation as an initial value. Thereby, for example, image reconstruction focusing on blood vessels in the brain can be performed. This is particularly effective for MRI image reconstruction of tubular objects such as blood vessels that are not filled inside.

According to the first solution of the present invention,
The processing unit stores a two-dimensional measurement image signal in which information on each cross section of the subject is imaged in measurement by magnetic resonance imaging and a predetermined parameter for the measurement, a storage unit, a measurement device, or an input unit Step to input from,
The processing unit obtains a first in-focus image having a focal plane with a desired depth position coordinate in the depth direction for the two-dimensional measurement image signal based on the input measurement image signal, and the obtained first Obtaining an absolute value amplitude image of one in-focus image;
The processing unit obtains a point spread function indicating a blur effect in which a video component other than the focal plane is superimposed on the focused image based on the input predetermined parameter, and an absolute value of the obtained point spread function A step of seeking an image;
The processing unit performs a Fourier transform on the obtained absolute value image of the obtained point spread function and the obtained absolute value amplitude image of the first focused image, obtains a ratio, and performs inverse Fourier transform to perform an inverse Fourier transform on the other than the focal plane. Obtaining an incoherent imaging approximate image from which a part or all of the video component is removed;
The processing unit convolves and integrates the obtained incoherent imaging approximate image and the point spread function to obtain a second focused image, and a ratio between the first focused image and the second focused image; Correlating the point spread function to obtain a correction function, correcting the incoherent image approximation image according to the obtained correction function to obtain a reconstructed image from which video components other than the focal plane are further removed;
The processing unit stores the obtained reconstruction image and / or incoherent imaging approximate image in the storage unit, or displays the display image on the display unit or outputs the output image to the output unit. An image reconstruction method is provided.

According to the second solution of the present invention,
The processing unit stores a two-dimensional measurement image signal in which information on each cross section of the subject is imaged in measurement by magnetic resonance imaging and a predetermined parameter for the measurement, a storage unit, a measurement device, or an input unit Step to input from,
The processing unit obtains a first in-focus image having a focal plane with a desired depth position coordinate in the depth direction for the two-dimensional measurement image signal based on the input measurement image signal, and the obtained first Obtaining an absolute value amplitude image of one in-focus image;
The processing unit obtains a point spread function indicating a blur effect in which a video component other than the focal plane is superimposed on the focused image based on the input predetermined parameter, and an absolute value of the obtained point spread function A step of seeking an image;
The processing unit performs a Fourier transform on the obtained absolute value image of the obtained point spread function and the obtained absolute value amplitude image of the first focused image, obtains a ratio, and performs inverse Fourier transform to perform an inverse Fourier transform on the other than the focal plane. Obtaining an incoherent imaging approximate image from which a part or all of the video component is removed;
The processing unit convolves and integrates the obtained incoherent imaging approximate image and the point spread function to obtain a second focused image, and a ratio between the first focused image and the second focused image; Correlating the point spread function to obtain a correction function, correcting the incoherent image approximation image according to the obtained correction function to obtain a reconstructed image from which video components other than the focal plane are further removed;
The processing unit causes the computer to execute a step of storing the obtained reconstructed image and / or incoherent imaging approximate image in the storage unit, displaying the image on the display unit, or outputting the output image to the output unit. An image reconstruction program is provided.

According to the third solution of the present invention,
A storage unit for storing a two-dimensional measurement image signal in which information of each cross section of a subject is imaged in a measurement by magnetic resonance imaging, and a predetermined parameter for the measurement;
A display unit and / or an output unit for displaying or outputting the determined predetermined image or information;
A processing unit for removing video components other than the focal plane from the first focused image based on the input measurement image signal,
The processing unit inputs a measurement image signal and a predetermined parameter for the measurement from the storage unit,
The processing unit obtains a first in-focus image with a desired depth position coordinate in the depth direction as a focal plane for the two-dimensional measurement image signal based on the input measurement image signal, and is obtained Obtaining an absolute amplitude image of the first focused image;
The processing unit obtains a point spread function indicating a blur effect in which video components other than the focal plane are superimposed on the focused image based on the input predetermined parameters, and the absolute value of the obtained point spread function Find the value image,
The processing unit performs a Fourier transform on the absolute value image of the obtained point spread function and the absolute value amplitude image of the obtained first in-focus image to obtain a ratio, and performs inverse filtering to perform an inverse Fourier transform, thereby performing a focal plane. Find an incoherent imaging approximate image from which part or all of the video components other than
The processing unit convolves and integrates the obtained incoherent image formation approximate image and the point spread function to obtain a second focused image, and calculates a ratio between the first focused image and the second focused image. Then, a correction function is obtained by performing a correlation operation with the point spread function, and an incoherent image approximation image is corrected according to the obtained correction function to obtain a reconstructed image from which video components other than the focal plane are further removed. ,
The processing unit stores the obtained reconstructed image and / or incoherent image approximation image in the storage unit, or displays the image on the display unit or outputs to the output unit. Provided.

  According to the present invention, there is provided an image reconstruction method and an image reconstruction program for removing a background image component superimposed on a focal plane image from a focal plane image group obtained in the Fresnel diffraction type imaging method (for example, see Non-Patent Document 1). Can be provided. Further, according to the present invention, by removing image information components other than the focal plane, for example, it is possible to easily identify a normal part and an abnormal part in image diagnosis.

  According to the present invention, MRI measurement time can be greatly shortened. According to the present invention, it is possible to shorten the MRI imaging time, which takes a long time, by performing three-dimensional information by two-dimensional data collection scanning. In addition, according to the present invention, it is possible to improve the throughput due to the increase in the number of MRI examinees per day, reduce the cost of consultation, and alleviate the mental distress of the subject in an obstructive environment. In addition, according to the present invention, it is possible to reduce the physiological effects (neural stimulation due to eddy current generation in the living body, increase in living body temperature) of the subject, and the utility is extremely large.

1. Outline The signal description formula of a signal (NMR signal) in the Fresnel diffraction type imaging method is an expression that approximates a three-dimensional diffraction wavefront distribution, and is given by the following expression.

  Here, υ (x ′, y ′) is an NMR signal to be measured, K ′ is a constant, ρ (x, y, z) is a subject distribution function, and b is a quadratic function of z = 0 plane. The coefficient of the changing magnetic field, τ is the application time of the magnetic field, α is the rate of change in the z direction of the magnetic field, (x ′, y ′) is the center coordinate of the quadratic function-like magnetic field, and γ is a constant specific to the nucleus (nucleus Magnetic rotation ratio).

  The signal of equation (1) is a two-dimensional signal acquisition, but information in the depth direction is also encoded in the signal by the z term (1 + αz) present on the shoulder of the exponential term. The Fresnel integral of Expression (1) is obtained by focusing on an arbitrary surface (z ′) in the depth direction by a back propagation process (inverse filtering) expressed by the following expression (first focused image) ρ ′ (x , Y, z) can be reconstructed.

_{F xy} and _F ^{xy -1} in Equation (2) is, x, Fourier transform and inverse Fourier transform about the _{y, k x,} k y represents a coordinate in the Fourier space. By setting a desired position coordinate z ′ (z = z ′) as the depth position coordinate z used in the calculation reconstruction, an image focused on an arbitrary z ′ plane can be obtained. That is, if the image is reconstructed by sequentially changing the focal plane in the depth direction from the two-dimensionally collected signal, the distribution of the subject can be known three-dimensionally. This is just like the process of adjusting the focus by changing the focal length of the lens when shooting with the camera. Since the imaging is two-dimensional, it is possible to know the three-dimensional distribution of the subject in a very short time.

  However, the image obtained in this case is clear in the focal plane image, but image components other than the focal plane are also superimposed on the image in a blurred form according to the distance from the focal plane. This is good for obtaining a sense of perspective, but there are cases in which it is not possible to distinguish between an image on the focal plane and an image other than the focal plane. The present embodiment relates to a method of removing an out-of-focus plane image component from a focal plane image in such a Fresnel diffraction type imaging method.

The focused image of Equation (2) can be approximately expressed by convolution integral as follows.

Here, β is a constant, and p _d (x, y, z) is a function that determines the blur effect (diffraction effect) when an image other than the focal plane is superimposed on the focal plane image, and is a point image distribution. function (PSF: Point Spread function) is a corresponding function (hereinafter, referred to as the point spread function _{p d).} The removal of the out-of-focus image is to obtain the subject distribution function ρ (x, y, z) from the equation (3). The subject distribution function at the desired position coordinate z ′ indicates a cross-sectional image of the subject at the position coordinate.

FIG. 1 is an explanatory diagram of general inverse filtering by complex number arithmetic. In the case of obtaining ρ (x, y, z) (FIG. 1A) from the convolution integral such as Expression (3), a method of performing image reconstruction by inverse filtering is generally known. That is, as shown in the following equation, the frequency spectrum (corresponding to FIG. 1C) of the focused image ρ ′ (x, y, z) is represented by the frequency spectrum of the point spread function p _d (x, y, z) ( This is a method of performing an inverse Fourier transform after dividing in FIG.
F ⁻¹ [F [ρ ′ (x, y, z)] / F [p _d (x, y, z)]]

As shown in FIG. 1B, the point spread function p _d (x, y, z) in the equation (4) is not a function distributed over the entire frequency space, and is a cutoff band in which no information exists in the frequency space. Exists. Information cannot originally be restored from where there is no information, and it is difficult to restore the spectrum of the subject distribution function ρ (x, y, z) by general inverse filtering, and the out-of-focus video component is removed. It is difficult. For example, it is difficult to restore the spectrum (broken line) to be obtained from the spectrum (solid line) in FIG.

  Therefore, in the present embodiment, the focal plane image component is extracted mainly using two methods in combination. First, the first method is inverse filtering of incoherent imaging approximation. The process of calculating a focused image by the Fresnel diffraction type imaging method is optically equivalent to coherent imaging in that the phase information of the wavefront is used. In this imaging, the focal plane image component cannot be obtained by the general inverse filtering as described above.

  FIG. 2 is an explanatory diagram of inverse filtering (inverse coherent imaging approximation inverse filtering) using an absolute value image according to the present embodiment.

As shown in FIG. 2, a case is considered where the phase information of the focused image is discarded and an absolute value amplitude image calculated from the amplitudes of the real part and the imaginary part is obtained. The absolute value amplitude image can be obtained from the square root of the sum of squares of the real part and the imaginary part of ρ ′ (x, y, z) represented by Expression (2).
√ ((real part of ρ ′) ² + (imaginary part of ρ ′) ² )
Note that although the real part and the imaginary part are not shown separately in Equation (2), the real part and the imaginary part are obtained because they are functions of complex numbers.

Absolute value amplitude images focused image is the convolution of the object function [rho (FIG. 2 (a)) and point spread absolute value image _p abs distribution function _{p d (x, y, z} ) ( FIG. 2 (b)) Can approximate an integral.
| Ρ ′ (x, y, z) | ≈ 1 / K ′ (ρ (x, y, z) * | p _d (x, y, z) |)
For example, when the subject has a small overlap of images on the plane and in the depth direction like a blood vessel, the absolute value image of the focused image approximates the convolution integral. As shown in FIG. 2B, the Fourier transform function of the absolute value image p _abs (x, y, z) of this point spread function is the Fourier transform of p _d (x, y, z) (FIG. 1 ( Unlike b)), the function is widely distributed in the frequency space. From this, the spectrum of the subject distribution function ρ can be restored by inverse filtering within the range where the approximation is established, and the out-of-focus video component can be removed. This inverse filtering by the absolute value image conversion of the focused image is similar to incoherent imaging in which only the intensity distribution of the wavefront is measured. Therefore, in the present embodiment, this is referred to as inverse filtering of incoherent imaging approximation.

The inverse filtering of this incoherent imaging approximation can be expressed by the following equation.
Note that ρ _inc obtained here is referred to as an incoherent imaging approximate image in the present embodiment.

  Next, as a second method, an image is estimated by a maximum likelihood method (Maximum Likelihood Method). Note that the iterative formula of the maximum likelihood method shown below is effective for image estimation (removal of out-of-focus components), and an image can be estimated only by the second method, but after the first method described above, The second method is more effective. For example, the calculation time of the maximum likelihood method can be shortened, and the calculation time can be shortened as a whole.

  Expression (6) is an iterative expression of deconvolution integral based on Bayesian theory (see Non-Patent Document 2). Image correction is performed by the effect of maximizing the likelihood in each step, and the original image component is estimated asymptotically.

For example, as a calculation procedure of the image restoration processing, first, an incoherent imaging approximate image ρ _inc obtained by performing inverse filtering of incoherent imaging approximation according to Expression (5) as an initial image is expressed by Expression (6). Substitute into ρ ^old (first image). Next, a focused image (second focused image) ρ ^old * p _d (x, y, z) obtained from the first image ρ ^old is calculated. Then, the obtained second in-focus image is compared with the first in-focus image ρ ′ (x, y, z) (true value) obtained according to the equation (2) based on the measurement image signal. , Correct the image. Since there is an error in the first image, the focused image ρ _old * p _d (x, y, z) calculated from the first image and the focused image ρ ′ (x, There is an error between y and z). For example, an error between the in-focus image ρ _old * p _d (x, y, z) calculated from the first image and the in-focus image ρ ′ (x, y, z) obtained from the measurement image signal is also included. Then, a correction function is calculated (calculation in parentheses in Equation (6)), and this is multiplied by ρ ^old to obtain ρ ^new (second image) (image correction processing). When the first image ρ ^old is multiplied by the correction function, the error of the first image ρ ^old is reduced. This multiplication is processing for correcting the image.

Thereafter, the obtained second image ρ ^new is substituted into the first image ρ ^old, and the processing after the calculation of the focused image ρ ^old * p _d (x, y, z) is repeated. For example, it can be repeated until the difference between ρ ^old and ρ ^new converges below a predetermined threshold, or can be repeated a predetermined number of times.

(Acquisition of measurement image by Fresnel diffraction imaging method)
Here, the outline of the measurement of the NMR signal by the Fresnel diffraction type imaging method will be described. As the measurement image signal (NMR signal) in the present embodiment, for example, an MRI image using a three-dimensional diffraction distribution can be used. For example, based on the NMR Fresnel diffraction type imaging method described in Non-Patent Document 1, an NMR signal obtained as follows can be used. This measurement is performed by an appropriate measurement apparatus such as an MRI apparatus, and a measurement signal can be input via an interface or the like.

FIG. 10 shows a pulse sequence for imaging by NMR Fresnel diffraction imaging. A quadratic function-like magnetic field, a quadratic function-like magnetic field scanning magnetic field, a sweep magnetic field, and a gradient magnetic field are changed as follows.
(A) Quadratic functional magnetic field: Δb (x, y, z) = b (1 + αz) (x ² + y ² )
(B) Quadratic functional magnetic field scanning magnetic field: Δb _0x (z) = (√ (1 + αz)) b _ox
(C) Sweep magnetic field: Δb _0y (z) = (1 + αz) b _oy
(D) Gradient magnetic field: ΔG _yz (y, z) = (1 + αz) × G _y y

After excitation of the spin system by the excitation pulse, the quadratic function-like magnetic field of (a) and the quadratic function-like magnetic field scanning magnetic field (hereinafter abbreviated as scanning magnetic field) of (b) are used in the phase encoding direction (for example, the x direction) ) For a time τ. The scanning magnetic field is a magnetic field used to scan the center position of the quadratic function-like magnetic field. When the scanning amount is x ′, the quadratic function-like magnetic field after scanning is
Δb = b (1 + αz) ((x′−x) ² + y ² )
And x ′ move the center position. Next, the magnetic field of (c) is applied by a sweep pulse, and the gradient magnetic field of (d) is applied by an inversion pulse in synchronization with the sweep magnetic field. The echo signal obtained at this time is a signal represented by Expression (1).

2. Hardware Configuration FIG. 3 is a hardware configuration diagram according to the present embodiment.
This hardware includes a processing unit 1, which is a central processing unit (CPU), an input unit 2, an output unit 3, a display unit 4, and a storage unit 5. Furthermore, you may have the interface 6 for connecting with a measuring apparatus, for example. The processing unit 1, the input unit 2, the output unit 3, the display unit 4, the storage unit 5, and the interface 6 are connected by appropriate connection means such as a star or a bus.

  The storage unit 5 includes a parameter file 51, a measurement image (NMR signal) file 52, a focused image file 53, a point spread function file 54, an incoherent imaging approximate image file 55, and a reconstructed image file 56. The storage unit 5 can also include a three-dimensional reconstructed image file 57. The above file may be omitted as appropriate.

  The parameter file 51 includes, for example, constants K ′ and β, a magnetic field coefficient b in which the intensity changes in a quadratic function shape of z = 0 plane, a magnetic field application time τ, a magnetic field z-direction change rate α, and a nucleus specific. A constant (nuclear magnetic rotation ratio) γ and the like are stored. These parameters can be stored in advance. In addition, parameters (for example, τ) for measurement of NMR signals may be appropriately input from an input unit or from a measurement device via an interface and stored in the parameter file 51.

The measurement image file 52 stores the NMR signal ν input from the measurement device or the like via the interface or input from the input unit. Focus image file 53, the point spread function file 54, respectively, focused image obtained based on the NMR signal [rho ', the point spread function p _d are stored. The incoherent imaging approximate image file 55 stores an incoherent imaging approximate image ρ _inc obtained by inverse filtering of incoherent imaging approximation according to the present embodiment.

The reconstructed image file 56 stores the reconstructed image ρ ^new from which the blur has been removed. The three-dimensional reconstructed image file 57 stores a three-dimensional image composed of a two-dimensional reconstructed image group.

3. Flow chart (image reconstruction)
FIG. 4 is a flowchart (1) of the present embodiment including inverse filtering of incoherent imaging approximation. FIG. 5 is a flowchart (2) of the present embodiment including the maximum likelihood calculation algorithm.

  First, the processing unit 1 inputs an NMR signal (measurement image signal) υ (x ′, y ′) and predetermined parameters from the measurement image file 52 and the parameter file 51 (S100). The processing unit 1 may input an NMR signal from an appropriate measuring device via the interface 6. The input parameters are, for example, constants K ′ and β, a magnetic field coefficient b whose intensity changes in a quadratic function of z = 0 plane, a magnetic field application time τ, a magnetic field z-direction change rate α, and a characteristic of the nucleus. Constant (nuclear magnetorotation ratio) γ and the like.

  Further, the processing unit 1 sets the depth position coordinate z to a desired coordinate z ′ (S101). For example, the processing unit 1 may read a predetermined coordinate z ′ from the parameter file 51 or may input it from the input unit 2. Further, the depth position coordinate z may be set according to a predetermined initial value and step size. A plurality of reconstructed images can be obtained by sequentially changing the depth position coordinate z, which will be described later.

  The processing unit 1 obtains a focused image ρ ′ focused on the set depth position coordinate z based on the input NMR signal (S102). For example, the processing unit 1 can obtain the focused image ρ ′ according to the above equation (2). The processing unit 1 may store the obtained focused image in the focused image file 53 in correspondence with the depth position coordinates.

The processing unit 1, based on the set depth position coordinate z and the input parameters, determining the point spread function p _d (S103). For example, the processing unit 1 can obtain a point spread function according to the above equation (4). The processing unit 1 may store the obtained point spread function in the point spread function file 54. Next, the processing unit 1 obtains an absolute value image p _abs of the obtained point spread function (S104).

The processing unit 1 assumes that the convolution integral between the subject distribution function ρ and the absolute value image p _abs of the obtained point spread function is the absolute value of the focused image ρ ′, and the absolute value of the obtained point spread function. By performing inverse filtering based on the absolute value amplitude image of the value image p _abs and the obtained focused image ρ ′, an incoherent imaging approximate image ρ _inc corresponding to the subject distribution function ρ is obtained (the above-mentioned first image Equivalent to the method). For example, the processing unit 1 can obtain an incoherent image formation approximate image according to the above-described equation (4).

As a specific procedure, first, the processing unit 1 obtains an absolute value amplitude image of the obtained focused image ρ ′, and performs a Fourier transform on the obtained absolute value amplitude image (S105). For example, the processing unit 1 can obtain an absolute value amplitude image by calculating the square root of the square sum of the real part and the imaginary part of the focused image ρ ′. Further, the processing unit 1 performs a Fourier transform on the absolute value image p _abs of the point spread function (S107).

Next, the processing unit 1 divides the absolute value amplitude image Fourier-transformed in Step S105 by the absolute value image of the point spread function Fourier-transformed in Step S107 (S109). Further, the processing unit 1 obtains an incoherent imaging approximate image ρ _inc by performing an inverse Fourier transform on the result divided in step S109 (S111). In the incoherent imaging approximate image ρ _inc obtained here, part or all of the video components other than the focal plane are removed.

Note that this embodiment, for example, in terms of using the absolute value image p _abs rather than using the point spread function p _d, performs inverse filtering the light intensity image having no phase information It is a technique.

The processing unit 1 stores the obtained incoherent image formation approximate image ρ _inc in the incoherent image formation approximate image file 55, or displays it on the display unit 4, outputs it to the output unit 3 (S113). The processing unit 1 may store an incoherent imaging approximate image corresponding to the set depth position coordinate z. Note that the process of step S113 may be omitted. In addition, although it is possible to omit the processes of steps S103 to S113, the convergence of iterative processes to be executed later can be accelerated by executing these steps.

Moving to FIG. 5, the processing unit 1 sets the obtained incoherent imaging approximate image ρ _inc as the first image ρ ^old (S201). For example, when the inverse filtering of incoherent imaging approximation is not performed (when steps S103 to S113 are omitted), the processing unit 1 uses the focused image ρ ′ obtained in step S102 as the first image ρ. ^old .

Processing unit 1, based on the first image [rho ^old, and the point spread function point spread function stored in the file 54 p _d, stored in the focused image file 53 focused image [rho 'and, A second image ρ ^new is obtained by the maximum likelihood method. For example, the processing unit 1 can obtain the second image ρ ^new according to the above equation (6).

As a specific procedure, first processing unit 1 reads out the point spread function _{p d} from the point spread function file 54, the first image [rho ^old and point spread function _{p d} convolution integral (S203). The processing unit 1, instead of reading the point spread function p _d, may be calculated again. Next, the processing unit 1 reads the focused image ρ ′ from the focused image file 53, and divides the focused image ρ ′ by the integration result (second focused image) in step S203 (S205). Note that the processing unit 1 may calculate again based on the NMR signal instead of reading the focused image ρ ′. Processor 1 obtains the division result obtained in step S205, the correction function by correlation calculation between the point spread function p _d (S207). Further, the processing unit 1 obtains a second image ρ ^new by multiplying the first image ρ ^old by the correction function obtained in step S207 (S209).

Next, the processing unit 1 determines convergence (S213). For example, the processing unit 1 obtains a difference between the obtained second image ρ ^new and the first image ρ ^old and compares the obtained difference with a predetermined value to determine convergence. . For example, the processing unit 1 can determine that the difference has converged when the obtained difference is equal to or less than a predetermined value (or smaller than that). In addition, instead of determining convergence, the processing unit 1 may determine whether the calculation has been performed a predetermined number of times (for example, 20 times) and determine whether to end the calculation.

When determining that the image has not converged (or does not end the calculation) (S213), the processing unit 1 sets the obtained second image ρ ^new as the first image ρ ^old (S215). Further, the processing unit 1 repeats the processes after step S203 described above based on the set new first image ρ ^old .

On the other hand, when it determines with having converged (or complete | finishing calculation) (S213), the process part 1 memorize | stores the calculated | required 2nd image (rho) ^new in the reconstruction image file 56, and / or The data is displayed on the display unit 4 and output to the output unit 3 (S217). The image obtained here is a reconstructed image from which image components other than the focal plane have been removed. Note that the processing unit 1 may store the reconstructed image in the reconstructed image file 56 in accordance with the set depth position coordinates. The processing unit 1 may store the obtained intermediate results and final results in the storage unit 5 as appropriate, or may read them out as appropriate.

(Multiple image reconstruction with varying depth position coordinates)
FIG. 6 is a flowchart for obtaining a plurality of reconstructed images. First, the processing unit 1 inputs an NMR signal (measurement image signal) and predetermined parameters from the parameter file 51 and the measurement image file 52 as described above (S100). Next, the processing unit 1 sets the depth position coordinate z (S301). For example, the processing unit 1 can set one of a plurality of coordinates stored in advance in the parameter file 51 as the depth position coordinate z. Further, the depth position coordinate z may be sequentially increased or decreased by a step size based on a predetermined initial value and a step size. Regarding the setting of the depth position coordinate z, one of a plurality of coordinates can be set by an appropriate method other than this.

The processing unit 1 executes the process of step S102 as described above. Further, the processing unit 1 obtains the point spread function _pd and its absolute value image p _abs (S302). Note that the processing in step S302 is the same as that in steps S103 and S104 described above, and a detailed description thereof will be omitted.

Next, the processing unit 1 obtains an incoherent imaging approximate image by inverse filtering processing (S303). The details of step S303 are the same as steps S105 to S111 described above. Next, the processing unit 1 obtains a reconstructed image ρ ^new by the maximum likelihood method (S305). Note that step S305 is the same as steps S201 to S215 described above, and a description thereof will be omitted. In addition, the processing unit 1 stores the obtained reconstructed image ρ ^new in the reconstructed image file 56 corresponding to the set depth position coordinate z (S307).

Next, the processing unit 1 determines whether a predetermined number of reconstructed images have been obtained (S309). In addition to determining whether a predetermined number of reconstructed images have been obtained, all of the plurality of coordinates stored in advance in the parameter file 51 have been processed, or the depth position coordinates are greater than or equal to a predetermined value (or less). It may be determined whether or not. In addition to this, it may be determined whether or not the process is to be terminated based on an appropriate criterion. If the processing unit 1 determines that a predetermined number of reconstructed images have not been obtained (S309), the processing unit 1 returns to step S301, and executes the processing after step S301 to obtain a ^new reconstructed image ρnew. In step S301, depth position coordinates different from the previous process are set. On the other hand, when the processing unit 1 determines that a predetermined number of reconstructed images have been obtained (S309), the processing ends. The processing unit 1 may display the focused image ρ ′ and the incoherent approximate image ρ _inc on the display unit 4 and output the output unit 3 at an appropriate timing.

(3D image playback)
FIG. 7 is a flowchart for reproducing a three-dimensional image based on a two-dimensional NMR signal. In addition, the thing of the same sign as the above-mentioned shows the process similar to the above-mentioned.

  First, the processing unit 1 executes the processes of steps S100, S301, S102, and S302 to S309 as described above. Note that the processing of each step is the same as described above, and thus detailed description thereof is omitted.

If the processing unit 1 determines that a predetermined number of reconstructed images have not been obtained (S309), the processing unit 1 returns to step S301 in the same manner as described above, and executes the subsequent processing. On the other hand, when the processing unit 1 determines that a predetermined number of reconstructed images have been obtained (S309), the processing unit 1 configures a three-dimensional image based on a plurality of reconstructed images stored in the reconstructed image file 56. (S311). Note that the processing unit 1 can form a three-dimensional image by an appropriate method. Further, the processing unit 1 stores the obtained three-dimensional image in the three-dimensional reconstructed image file 57 and / or displays it on the display unit 4, outputs it to the output unit 3 (S313). The processing unit 1 may display the focused image ρ ′ and the incoherent approximate image ρ _inc on the display unit 4 and output the output unit 3 at an appropriate timing.
Thus, in the embodiment, a three-dimensional image can be reconstructed based on one two-dimensional measurement image.

4). Simulation Result FIG. 8 shows the result of a numerical simulation using a wire frame. In this example, the three-dimensional subject model is composed of 128 × 128 × 8 pixels. For example, for the illustrated three-dimensional object model (imaging model), depth position information z is set so as to focus on each of the back surface, center surface, and front surface of the model, and in-focus image and incoherent image approximation FIG. 8 shows the result of calculating the image and the reconstructed image.

  FIG. 8A is an image (focused image) focused on each of the back, center, and front surfaces of the model. Although the focal plane image is clear, the images of the planes before and after the focal plane are blurred and superimposed.

  FIG. 8B is an inverse filtered image of incoherent imaging approximation. Although the out-of-focus image is removed to some extent, an error remains because the incoherent imaging approximation is not completely established.

  FIG. 8C is an image estimated image by the maximum likelihood method. It can be seen that the image is recovered well. The image shown in the figure is a result of repeated processing 20 times. In this figure, the model image intensity is constant, but a halftone image can be recovered similarly. Image filtering can be estimated without performing incoherent imaging approximation inverse filtering in most cases as the first method, but using incoherent imaging approximation inverse filtering together has the effect of speeding up convergence. , It is possible to shorten the iteration time.

  FIG. 9 shows a three-dimensional image reconstruction simulation result of a model simulating blood vessels (blood vessel simulation numerical model). FIG. 9B shows an image obtained by simulating the three-dimensional model of FIG. 9A and focusing on each cross section from the obtained signal. As shown in FIG. 9B, images out of the focal plane are blurred and superimposed. In FIG. 9, four cross-sectional images are shown for each of the in-focus image, the incoherent imaging approximate image, and the reconstructed image. However, in the simulation, 32 cross-sectional images are estimated and a three-dimensional image is obtained. Was reconfigured.

  FIG. 9C is an inverse filtering image of incoherent imaging approximation. The result of applying this inverse filtered image as the initial image to the maximum likelihood method is the image of FIG. 9D. Further, FIG. 9E shows the result of constructing a three-dimensional image from the cross-sectional image group obtained here. The three-dimensional image can be constructed by an appropriate method based on the obtained cross-sectional image group. It can be seen that the reconstructed three-dimensional blood vessel image (FIG. 9E) matches the image set in the model (FIG. 9A) with extremely high accuracy.

5). APPENDIX The image reconstruction method of the present invention includes an image reconstruction program for causing a computer to execute each procedure, a computer-readable recording medium recording the image reconstruction program, and an image reconstruction program in an internal memory of the computer. It can be provided by a loadable program product, a computer such as a server including the program, and the like.

  The present invention can be used in, for example, an industry for manufacturing an apparatus using a magnetic resonance imaging method using a Fresnel diffraction method, and other industries related to image sensing.

Explanatory drawing of general inverse filtering by complex number operation. Explanatory drawing of the inverse filtering (incoherent image formation approximate inverse filtering) by the absolute value image of this Embodiment. The hardware block diagram regarding this Embodiment. The flowchart (1) of this Embodiment including an incoherent image formation approximate inverse filtering. The flowchart (2) of this Embodiment containing the calculation algorithm of a maximum likelihood method. Flow chart for obtaining multiple reconstructed images The flowchart which reproduces | regenerates a 3-dimensional image based on a 2-dimensional NMR signal. Results of numerical simulation using wire frame. Simulation results of a model that simulates a blood vessel. Pulse sequence for imaging by NMR Fresnel diffraction imaging.

Explanation of symbols

DESCRIPTION OF SYMBOLS 1 Processing part 2 Input part 3 Output part 4 Display part 5 Memory | storage part 6 Interface 51 Parameter file 52 Measurement image (NMR signal) file 53 Focus image file 54 Point image distribution function file 55 Incoherent imaging approximate image file 56 Reconstruction Image file 57 3D reconstruction image file

Claims (12)

The processing unit stores a two-dimensional measurement image signal in which information on each cross section of the subject is imaged in measurement by magnetic resonance imaging and a predetermined parameter for the measurement, a storage unit, a measurement device, or an input unit Step to input from,
The processing unit obtains a first in-focus image having a focal plane with a desired depth position coordinate in the depth direction for the two-dimensional measurement image signal based on the input measurement image signal, and the obtained first Obtaining an absolute value amplitude image of one in-focus image;
The processing unit obtains a point spread function indicating a blur effect in which a video component other than the focal plane is superimposed on the focused image based on the input predetermined parameter, and an absolute value of the obtained point spread function A step of seeking an image;
The processing unit performs a Fourier transform on the obtained absolute value image of the obtained point spread function and the obtained absolute value amplitude image of the first focused image, obtains a ratio, and performs inverse Fourier transform to perform an inverse Fourier transform on the other than the focal plane. Obtaining an incoherent imaging approximate image from which a part or all of the video component is removed;
The processing unit convolves and integrates the obtained incoherent imaging approximate image and the point spread function to obtain a second focused image, and a ratio between the first focused image and the second focused image; Correlating the point spread function to obtain a correction function, correcting the incoherent image approximation image according to the obtained correction function to obtain a reconstructed image from which video components other than the focal plane are further removed;
The processing unit stores the obtained reconstructed image and / or incoherent imaging approximate image in the storage unit, or displays the image on the display unit or outputs the image to the output unit. The processing unit sequentially sets a plurality of depth position coordinates;
The image reconstructing method according to claim 1, further comprising: a step of repeatedly obtaining a reconstructed image for each of the set depth position coordinates. The processing unit sequentially sets a plurality of depth position coordinates;
A processing unit that repeats obtaining a reconstructed image for each of the set depth position coordinates; and
A processing unit configured to form a three-dimensional image corresponding to the subject based on the plurality of reconstructed images obtained;
The image reconstructing method according to claim 1, wherein the processing unit further includes a step of storing the configured three-dimensional image in the storage unit, displaying on the display unit, or outputting to the output unit. The step of obtaining the incoherent imaging approximate image includes:
The processing unit performs a Fourier transform on the obtained absolute value amplitude image of the first focused image;
The processing unit performs a Fourier transform on the absolute value image of the obtained point spread function,
The processing unit divides the absolute value amplitude image of the converted first focused image by the absolute value image of the converted point spread function;
The image processing method according to claim 1, wherein the processing unit includes a step of performing an inverse Fourier transform on the division result to obtain an incoherent imaging approximate image. The image reconstruction method according to claim 1, wherein in the step of obtaining the incoherent imaged approximate image, the processing unit obtains an incoherent imaged approximate image according to the following equation.
P _d : Point spread function, P _abs : Absolute value image of point spread function, ρ ′: First focused image, ρ _inc : Incoherent imaging approximate image, F: Fourier transform, F ⁻¹ : Inverse Fourier transform, z: depth position coordinate The image reconstruction method according to claim 1, wherein the step of obtaining an absolute value image of the point spread function includes a processing unit obtaining a point spread function according to the following equation.
Here, b: coefficient of the magnetic field whose intensity changes in a quadratic function on the z = 0 plane, τ: application time of the magnetic field, γ: a constant specific to the nucleus (nuclear magnetorotation ratio), β: constant, z : Depth position coordinates The step of obtaining the reconstructed image includes:
The processing unit sets the obtained incoherent imaging approximate image as a first image,
The processing unit convolves and integrates the first image and the obtained point spread function to obtain a second focused image;
The processing unit divides the first focused image by the second focused image;
The processing unit correlates the division result of the dividing step and the obtained point spread function to obtain a correction function;
The processing unit multiplies the first image and the correction function to obtain a reconstructed image;
The processing unit newly sets the obtained reconstructed image as the first image until the difference between the reconstructed image and the first image converges to a predetermined value or less or repeats a predetermined number of times. The method of reconstructing an image according to claim 1, further comprising the step of repeatedly obtaining a reconstructed image based on the first image. The step of obtaining the reconstructed image includes:
The processing unit sets the obtained incoherent imaging approximate image as a first image,
The processing unit obtains a reconstructed image according to the following formula based on the first image;
The processing unit newly sets the obtained reconstructed image as the first image until the difference between the reconstructed image and the first image converges to a predetermined value or less or repeats a predetermined number of times. The method of reconstructing an image according to claim 1, further comprising the step of repeatedly obtaining a reconstructed image based on the first image.
Ρ ′: first focused image, P _d : point spread function, ρ ^old : first image, ρ ^new : reconstructed image, z: depth position coordinates   2. The measurement image signal is an NMR signal, a measurement image signal obtained by a magnetic resonance imaging method using a Fresnel diffraction method, or a measurement image signal obtained by a magnetic resonance imaging method using a three-dimensional diffraction distribution. Image reconstruction method. 2. The image reconstruction method according to claim 1, wherein the step of obtaining the absolute value amplitude image of the first focused image includes a processing unit obtaining a first focused image according to the following equation.
Here, ρ ′ (x, y, z): first focused image, υ (x ′, y ′): measurement image signal, K ′: constant, b: quadratic function of z = 0 plane Coefficient of magnetic field whose intensity changes, τ: application time of the magnetic field, α: rate of change in the z direction of the magnetic field, (x ′, y ′): central coordinates of quadratic functional magnetic field, γ: constant specific to the nucleus ( Nuclear magnetic rotation ratio), z: depth position coordinates, k _x , k _y : coordinates in Fourier space, F _xy : Fourier transform for x and y, F _xy ⁻¹ : inverse Fourier transform for x and y. The processing unit stores a two-dimensional measurement image signal in which information on each cross section of the subject is imaged in measurement by magnetic resonance imaging and a predetermined parameter for the measurement, a storage unit, a measurement device, or an input unit Step to input from,
The processing unit obtains a first in-focus image having a focal plane with a desired depth position coordinate in the depth direction for the two-dimensional measurement image signal based on the input measurement image signal, and the obtained first Obtaining an absolute value amplitude image of one in-focus image;
The processing unit obtains a point spread function indicating a blur effect in which a video component other than the focal plane is superimposed on the focused image based on the input predetermined parameter, and an absolute value of the obtained point spread function A step of seeking an image;
The processing unit performs a Fourier transform on the obtained absolute value image of the obtained point spread function and the obtained absolute value amplitude image of the first focused image, obtains a ratio, and performs inverse Fourier transform to perform an inverse Fourier transform on the other than the focal plane. Obtaining an incoherent imaging approximate image from which a part or all of the video component is removed;
The processing unit convolves and integrates the obtained incoherent imaging approximate image and the point spread function to obtain a second focused image, and a ratio between the first focused image and the second focused image; Correlating the point spread function to obtain a correction function, correcting the incoherent image approximation image according to the obtained correction function to obtain a reconstructed image from which video components other than the focal plane are further removed;
The processing unit causes the computer to execute a step of storing the obtained reconstructed image and / or incoherent imaging approximate image in the storage unit, displaying the image on the display unit, or outputting the output image to the output unit. Image reconstruction program. A storage unit for storing a two-dimensional measurement image signal in which information of each cross section of a subject is imaged in a measurement by magnetic resonance imaging, and a predetermined parameter for the measurement;
A display unit and / or an output unit for displaying or outputting the determined predetermined image or information;
A processing unit for removing video components other than the focal plane from the first focused image based on the input measurement image signal,
The processing unit inputs a measurement image signal and a predetermined parameter for the measurement from the storage unit,
The processing unit obtains a first in-focus image with a desired depth position coordinate in the depth direction as a focal plane for the two-dimensional measurement image signal based on the input measurement image signal, and is obtained Obtaining an absolute amplitude image of the first focused image;
The processing unit obtains a point spread function indicating a blur effect in which video components other than the focal plane are superimposed on the focused image based on the input predetermined parameters, and the absolute value of the obtained point spread function Find the value image,
The processing unit performs a Fourier transform on the absolute value image of the obtained point spread function and the absolute value amplitude image of the obtained first in-focus image to obtain a ratio, and performs inverse filtering to perform an inverse Fourier transform, thereby performing a focal plane. Find an incoherent imaging approximate image from which part or all of the video components other than
The processing unit convolves and integrates the obtained incoherent image formation approximate image and the point spread function to obtain a second focused image, and calculates a ratio between the first focused image and the second focused image. Then, a correction function is obtained by performing a correlation operation with the point spread function, and an incoherent image approximation image is corrected according to the obtained correction function to obtain a reconstructed image from which video components other than the focal plane are further removed. ,
The processing unit stores the obtained reconstructed image and / or incoherent imaging approximate image in the storage unit, displays the image on the display unit, or outputs the image to the output unit.