KR102329938B1 - Method for processing conebeam computed tomography image using artificial neural network and apparatus therefor - Google Patents

Abstract

Disclosed are a method and apparatus for processing cone beam tomography images using a neural network. An image processing method according to an embodiment of the present invention includes receiving cone-beam computed tomography data; and reconstructing an image of the cone-beam computed tomography data using a neural network that removes cone-beam artifacts, wherein reconstructing the image includes a differential backprojection (DBP) algorithm. obtaining initial images in at least two directions from the cone beam computed tomography data using an analytic algorithm; and reconstructing an image of the cone-beam computed tomography data from the initial images using the neural network, wherein the reconstructing of the image includes a final image of the reconstructed image using a spectral blending technique. A restored image can be created.

Description

Cone beam tomography image processing method and apparatus using neural network

The present invention relates to a method and apparatus for processing a cone-beam tomography image, and more particularly, to an image processing method capable of reconstructing a cone-beam tomography image into a high-quality image using a neural network that removes cone-beam artifacts and to the device.

Cone-Beam Computed Tomography (CBCT) with a circular orbit is one of the frequently used CT systems due to its simple geometry. Since such CBCT can acquire a high-resolution projection image with a relatively simple scanner shape, it is often used for interventional imaging and dental CT. In general, a reconstruction algorithm called Feldkamp, Davis and Kress (FDK) is generally used to reconstruct an image using projection data obtained from a CBCT system. However, in the case of the FDK algorithm, as the cone-angle of the CBCT system increases, image noise called a cone-beam artifact becomes severe, which causes distortion of the shape of the reconstructed image.

Mathematically, cone beam artifacts are caused by defects inherent in circular trajectories, which do not satisfy Tuy's condition. This appears as a missing spectral component at a specific frequency determined by the scanner geometry in the Fourier domain. Therefore, accurate removal of cone beam artifacts may not be possible without additional prior information.

In order to remove such image distortion or image noise, various model-based iterative reconstruction algorithms (Model-Based Iterative Reconstruction, MBIR) have been generally used, but the CT projection operation and back-projection operation, which require a lot of computation, are repeatedly performed. There are limits to what must be done.

Recently, studies on computed tomography reconstruction techniques using neural networks have begun to be actively conducted. It has already been verified that low-dose CT, sparse-view CT, and internal tomography CT show a high level of restoration performance using neural networks. In the case of a restoration technique using a neural network, it is known that it shows excellent performance despite the very short restoration execution time. In addition, recent theoretical studies have shown that deep convolutional neural networks with encoder-decoder structures are related to new frame extensions using combinatorial convolutional frames.

Embodiments of the present invention provide an image processing method and apparatus capable of reconstructing a cone-beam tomography image into a high-quality image using a neural network that removes cone-beam artifacts.

Embodiments of the present invention restore a cone-beam tomography image to a high-quality image by removing residual noise using a spectral blending technique in which a cone-beam image from which cone-beam artifacts are removed by a neural network is removed as an input. Provided are an image processing method and an apparatus for the same.

An image processing method according to an embodiment of the present invention includes receiving cone-beam computed tomography data; and reconstructing an image of the cone-beam computed tomography data using a neural network that removes cone-beam artifacts.

The reconstructing of the image may include: obtaining initial images in at least two directions from the cone-beam computed tomography data using an analytic algorithm including a differential backprojection (DBP) algorithm; and reconstructing an image of the cone-beam computed tomography data from the initial images by using the neural network.

The restoring of the image includes obtaining an initial image in a first direction and an initial image in a second direction orthogonal to the first direction using a differential inverse projection algorithm, and using the neural network to obtain a cone beam in the first direction. The first image from which the artifacts are removed and the second image from which the cone beam artifacts in the second direction are removed may be reconstructed.

The reconstructing of the image may include generating a final reconstructed image of the cone-beam computed tomography data by using a spectral blending technique in which the first image and the second image are input.

The reconstructing of the image may include generating a final reconstructed image of the reconstructed image by using a spectral blending technique.

The neural network may include a neural network based on a convolution framelet.

The neural network may include a multi-resolution neural network including a pooling layer and an unpooling layer.

The neural network may include a bypass connection from the pooling layer to the unpooling layer.

An image processing apparatus according to an embodiment of the present invention includes: a receiver configured to receive cone-beam computed tomography data; and a restoration unit for reconstructing an image of the cone-beam computed tomography data using a neural network that removes cone-beam artifacts.

The restoration unit acquires initial images in at least two directions from the cone-beam computed tomography data using an analytic algorithm including a differential backprojection (DBP) algorithm, and uses the neural network to obtain the initial images An image of the cone-beam computed tomography data may be reconstructed from the images.

The restoration unit obtains an initial image in a first direction and an initial image in a second direction orthogonal to the first direction by using a differential back-projection algorithm, and removes the cone beam artifact in the first direction using the neural network. The first image and the second image from which the cone beam artifact in the second direction is removed may be reconstructed.

The reconstructor may generate a final reconstructed image for the cone-beam computed tomography data by using a spectral blending technique in which the first image and the second image are input.

The reconstructor may generate a final reconstructed image of the reconstructed image by using a spectral blending technique.

The neural network may include a neural network based on a convolution framelet.

The neural network may include a multi-resolution neural network including a pooling layer and an unpooling layer.

The neural network may include a bypass connection from the pooling layer to the unpooling layer.

According to embodiments of the present invention, a cone beam artifact is removed with a cone beam tomography image using a neural network that removes the cone beam artifact, and a spectral blending technique in which the cone beam image from which the cone beam artifact has been removed is used as an input. By removing the noise, the cone-beam tomography image can be restored to a high-quality image.

Computed tomography (CT) is a technique for reconstructing a tomography image of an experimenter using a projection image obtained by projecting X-rays onto an experimenter. In particular, cone-beam computed tomography (CBCT) is used in various fields as 3D CBCT for dentistry, 3D CBCT for plastic surgery, and 3D CBCT image for monitoring the patient's treatment progress because the track of the equipment used is simple. According to embodiments of the present invention, by performing deterioration correction through a differential back-projection (DBP) image, the performance may show stable performance regardless of uses such as dentistry or plastic surgery. Therefore, the present invention can be applied to various 3D CBCT systems.

3D CBCT for dentistry and 3D CBCT for plastic surgery have a high frequency of use, and a lot of equipment has already been introduced. As such, in a situation in which the infrastructure for 3D CBCT is sufficiently equipped, a core application such as the present invention has not been developed. Therefore, the present invention can be introduced immediately because it is sufficiently established in the infrastructure that can be directly applied to 3D CBCT for dental and plastic surgery, and through this, it is possible to provide high-quality images.

1 is a flowchart illustrating an image processing method according to an embodiment of the present invention.
2 shows an exemplary diagram for differential inverse-projection domain deep learning for explaining the method of the present invention.
3 shows an exemplary view of a coordinate system of a factorized cone beam geometry.
4 is a diagram illustrating an example of a source geometry and a filtering direction on a plane of interest.
5 is a diagram illustrating an example for explaining a missing frequency region.
6 is a diagram illustrating an example for explaining a spectral blending technique.
7 is a diagram illustrating an example of a disk phantom.
8 is a diagram illustrating an example of input and output of a neural network that removes cone beam artifacts.
9 shows an exemplary diagram of a deep neural network architecture in the method of the present invention.
10 shows an exemplary diagram for convergence of the objective function according to the number of epochs.
11 and 12 show exemplary views comparing restoration performance by an iterative restoration technique (total variation), FDK-CNN, and the method of the present invention.
13 illustrates a configuration of an image processing apparatus according to an embodiment of the present invention.

Advantages and features of the present invention and methods of achieving them will become apparent with reference to the embodiments described below in detail in conjunction with the accompanying drawings. However, the present invention is not limited to the embodiments disclosed below, but will be implemented in various different forms, and only these embodiments allow the disclosure of the present invention to be complete, and common knowledge in the art to which the present invention pertains It is provided to fully inform those who have the scope of the invention, and the present invention is only defined by the scope of the claims.

The terminology used herein is for the purpose of describing the embodiments, and is not intended to limit the present invention. As used herein, the singular also includes the plural unless specifically stated otherwise in the phrase. As used herein, "comprises" and/or "comprising" refers to the presence of one or more other components, steps, operations and/or elements mentioned. or addition is not excluded.

Unless otherwise defined, all terms (including technical and scientific terms) used herein may be used with the meaning commonly understood by those of ordinary skill in the art to which the present invention belongs. In addition, terms defined in a commonly used dictionary are not to be interpreted ideally or excessively unless specifically defined explicitly.

Hereinafter, preferred embodiments of the present invention will be described in more detail with reference to the accompanying drawings. The same reference numerals are used for the same components in the drawings, and repeated descriptions of the same components are omitted.

Embodiments of the present invention have a gist of restoring a cone-beam tomography image to a high-quality image by removing the cone-beam artifact with the cone-beam tomography image using a neural network that removes the cone-beam artifact.

Here, the present invention can restore a cone beam tomography image to a high quality image by removing residual noise using a spectral blending technique that takes a cone beam image from which the cone beam artifact has been removed as an input.

In this case, the neural network used in the present invention may include a convolution framelet-based neural network, and includes a multi-resolution neural network including a pooling layer and an unpooling layer. can do. Furthermore, the multi-resolution neural network may include a bypass connection from the pooling layer to the unpooling layer.

The neural network of the present invention inputs initial images for at least two directions obtained from cone-beam computed tomography images using an analytic algorithm including a differential backprojection (DBP) algorithm for cone-beam computed tomography images. , and by removing cone-beam artifacts for the initial images, it is possible to reconstruct a high-quality image of the cone-beam computed tomography image.

1 is a flowchart illustrating an image processing method according to an embodiment of the present invention, and FIG. 2 is a diagram illustrating an example of differential inverse-projection domain deep learning for explaining the method of the present invention.

1 and 2 , the image processing method according to an embodiment of the present invention includes the steps of receiving cone-beam computed tomography data (S110), and adding the cone-beam computed tomography data to the cone-beam computed tomography data using a neural network that removes cone-beam artifacts. It includes a step of reconstructing an image of the . (S120) and a step of generating a final reconstructed image of the reconstructed image using a spectral blending technique (S130).

Here, in step S120, the internal computed tomography data are first images in at least two directions from the cone beam computed tomography data using a preset analytic algorithm, for example, a differential backprojection (DBP) algorithm. For example, first images for each of the coronal and sagittal directions are acquired, and from the first images for at least two or more directions using a neural network that takes the first images (CNN input) as input. By outputting an image (CNN output) from which the cone-beam artifact is removed, an image for the cone-beam computed tomography data may be reconstructed. In this case, in step S120, initial images in two directions orthogonal to each other may be obtained, or initial images in two non-orthogonal directions may be obtained, and cone beam artifacts of each of the initial images are removed by using a neural network, so that the cone beam Reconstructed images in each of two directions from which artifacts are removed may be obtained. That is, as shown in FIG. 2 , a differential inverse projection operation is performed to generate projection data as a differential inverse projection (DBP) image in two different directions, for example, a coronal direction and a sagittal direction, and DBP images in both directions. A reconstructed image is generated by applying a neural network to each.

The neural network used in the present invention may include a convolution framelet-based neural network, and may include a multi-resolution neural network including a pooling layer and an unpooling layer. .

Here, the convolutional framelet may refer to a method of expressing an input signal using a local basis and a non-local basis.

Furthermore, the neural network may include a bypass connection from the pooling layer to the unpooling layer.

When the cone beam artifacts of the original images in at least two directions are removed through the neural network in step S120, spectral blending is performed by inputting spectral blending inputs in the at least two directions in which the cone beam artifacts are removed. A final reconstructed image is generated using the technique (S130).

Here, in step S130, a 3D image is reconstructed through the initial images for each direction by using each of the initial images for at least two directions from which the cone beam artifact is removed, and a slice in a different direction for each of the reconstructed 3D images A final reconstructed image may be generated by removing residual noise using a spectral blending technique using an image pair as an input. That is, in step S130, the noise remaining in the slice image pair in the other direction is removed to obtain a slice image in the other direction from which the remaining noise is removed, for example, a slice image in the third direction, and the slice image in the third direction from which the noise is removed. High-quality 3D images of cone-beam computed tomography data can be generated using

The method according to the present invention will be described with reference to FIGS. 3 to 12 .

notation

의 점이라 하면 x, y 및 z는 환자에 대한 직교 좌표를 의미한다. 원형 궤적의 경우 엑스선 소스(source)가 물체 f(x)를 중심으로 회전하면 소스 궤적은 아래 <수학식 1>과 같이 나타낼 수 있다.x = (x, y, z)

For a point in x, y and z mean the Cartesian coordinates for the patient. In the case of a circular trajectory, when the X-ray source rotates around the object f(x), the source trajectory can be expressed as in Equation 1 below.

[Equation 1]

Here, R may mean a radius of a circular scan, and symbol ┳ may mean transpose of a matrix or vector.

_{An X-ray transformation D f} that maps f(x) to a set of line integrals may be defined as in Equation 2 below.

[Equation 2]

의 벡터를 의미할 수 있다.where θ is the unit sphere

can mean a vector of .

The 3D inverse Fourier transform can be expressed as in Equation 3 below.

[Equation 3]

는 f(x)의 푸리에 스펙트럼을 의미하고,

은 각주파수를 의미하며,

을 의미할 수 있다.here,

is the Fourier spectrum of f(x),

is the angular frequency,

can mean

difference back projection (DBP)

의 투영 데이터를 이용하는 점(point)

에서 차분 역투영(DBP)은 아래 <수학식 4>와 같이 나타낼 수 있다.For a given X-ray source trajectory a(λ), X-ray source trajectory a(λ),

point using projection data of

In , the differential back-projection (DBP) can be expressed as in Equation 4 below.

[Equation 4]

Here, 1/||x-a(λ)|| may mean magnification factor dependent weighting.

One of the most fundamental properties of DBP is its relationship to the Fourier transform. Specifically, the DBP data of Equation 4 may be expressed as in <Equation 5> to <Equation 7> below.

[Equation 5]

[Equation 6]

[Equation 7]

라는 추가적인 항이다. 이것은 DBP가 f(x)의 필터링된 버전을 제공한다는 것을 의미한다. 필터

의 스펙트럼 또한 x에 의존하기 때문에 해당 필터는 공간적으로 변화하는 필터이다.Compared with the Fourier formula of Equation 3, the main difference of Equation 5 is

is an additional term. This means that DBP provides a filtered version of f(x). filter

Since the spectrum of is also dependent on x, the filter is a spatially varying filter.

Corn Beam geometric printing Factorized Representation of Conebeam Geometry)

In the case of cone beam CT with a circular trajectory, a factorization method has been proposed in previous studies, and similar ideas have been proposed in other studies as well. The present invention can utilize the core idea of this factorization method.

가 관심 평면을 의미한다 가정한다. 인자화 방법의 주요 목표는 3D 복원 문제를 관심 평면에서 연속적인 2D 문제로 전환하는 것이다. 특히 관심 평면

에서 가상의 코드선(chord line)과 가상의 소스 위치를 정의한다. 가상 코드선 좌표계는 가상 코드선 방향 e, z축 e_z 및 수직축 e┴에 의해 정의된다. 가상의 소스 a_v(λ⁺)와 a_v(λ^-)는 아래 <수학식 8>과 같이 간단히 계산될 수 있다.3 shows an exemplary view of a coordinate system of a factorized cone beam geometry, and shows a geometric structure of a cone beam CT having a circular trajectory. Here, consider a plane parallel to the z-axis and intersecting the source trajectory at two positions ^{a(λ -} ) and a(λ ^{+ ).}

Assume that is the plane of interest. The main goal of the factorization method is to transform the 3D reconstruction problem into a continuous 2D problem in the plane of interest. plane of particular interest

defines a virtual chord line and a virtual source location. The virtual code line coordinate system is defined by the virtual code line direction e, the z-axis e _z and the vertical axis e┴. The virtual sources a _v (λ ⁺ ) and a _v (λ ^- ) can be simply calculated as shown in Equation 8 below.

[Equation 8]

Then, the Cartesian coordinate x is converted into a new coordinate (t, s, z) as shown in Equation 9 below.

[Equation 9]

If the notation is slightly abused, the object density and DBP data in the plane of interest P of fixed s can be expressed as a 2D function as shown in Equation 10 below.

[Equation 10]

의 첫 번째 항과 두 번째 항은 각각 α(λ^-, x)와 α(λ⁺, x) 방향에 따른 힐버트(Hilbert) 변환에 해당하기 때문에 기존 연구(F. Dennerlein, F. Noo, H. Schondube, G. Lauritsch, and J. Hornegger, "A factorization approach for cone-beam reconstruction on a circular short-scan," IEEE Transactions on Medical Imaging, vol. 27, no. 7, pp. 887-896, 2008.)를 통해 상기 수학식 5를 아래 <수학식 11>과 같이 나타낼 수 있다.of Equation 6 above

The first and second terms of α(λ ^- , x) and α(λ ⁺ , Because x) corresponds to the Hilbert transform along the direction, previous studies (F. Dennerlein, F. Noo, H. Schondube, G. Lauritsch, and J. Hornegger, “A factorization approach for cone-beam reconstruction on a circular Short-scan," IEEE Transactions on Medical Imaging, vol. 27, no. 7, pp. 887-896, 2008.) through Equation 5 can be expressed as <Equation 11> below.

[Equation 11]

과

선을 따른 좌표를 의미할 수 있다. 이는 DBP 도메인에서 디컨볼루션 문제를 해결함으로써 콘빔 재구성 문제를 해결할 수 있음을 의미한다.Here, h _H (t) means the Hilbert transform, and z ₁ (τ) and z ₂ (τ) are as shown in FIG. 4 .

class

It may mean coordinates along a line. This means that the cone beam reconstruction problem can be solved by solving the deconvolution problem in the DBP domain.

에서 디컨볼루션을 위한 인코더-디코더 CNN

at Encoder-decoder CNNs for deconvolution

의 g(t, z)에서 원본 영상 f(t, z)를 복구할 수 있다고 제안하지만, 두 가지 기술적 어려움이 있다. 첫째, t에 의존하는 z₁(t) 및 z₂(t) 좌표 때문에 알려지지 않은 영상의 그리드는 표준 직교 그리드가 아니다. 실제로 f(t, z₁(t))와 f(t, z₂(t))는 새로운 좌표계 상 f(t, z)의 변형된 영상(deformed image)으로 볼 수 있다. 따라서, 컨볼루션 필터는 공간적으로 다르기 때문에 표준 디컨볼루션 방법은 작동하지 않는다. 둘째, 각 DBP 데이터 g(t, z)는 f(t, z₁(t))와 f(t, z₂(t))라는 두 개의 변형된 영상으로부터 기여하기 때문에 디컨볼루션 문제는 대단히 부적절하게 정의 되어진다. 따라서, 기존 연구(F. Dennerlein, F. Noo, H. Schondube, G. Lauritsch, and J. Hornegger, "A factorization approach for cone-beam reconstruction on a circular short-scan," IEEE Transactions on Medical Imaging, vol. 27, no. 7, pp. 887-896, 2008.)는 알려지지 않은 영상에 대하여 원래의 고정된 직교 그리드에 대하여 통합 이산 후 정규화된 행렬 반전 접근법을 제안하였다.Equation (11) is the deconvolution algorithm for each plane of interest.

We propose that we can recover the original image f(t, z) from g(t, z), but there are two technical difficulties. First, the grid of the unknown image is not a standard orthogonal grid because of the t- _{dependent z 1} (t) and z _{2 (t) coordinates.} In fact, f(t, z ₁ (t)) and f(t, z ₂ (t)) can be viewed as deformed images of f(t, z) on the new coordinate system. Therefore, the standard deconvolution method does not work because the convolution filters are spatially different. Second, since each DBP data g(t, z) _{contributes from two transformed images f(t, z 1} (t)) and f(t, z ₂ (t)), the deconvolution problem is very inappropriate. is well defined Therefore, previous studies (F. Dennerlein, F. Noo, H. Schondube, G. Lauritsch, and J. Hornegger, “A factorization approach for cone-beam reconstruction on a circular short-scan,” IEEE Transactions on Medical Imaging, vol. (27, no. 7, pp. 887-896, 2008.) proposed a normalized matrix inversion approach after integration discrete on the original fixed orthogonal grid for unknown images.

In general, the deconvolution algorithm for Equation 11 can be expressed as a mapping as shown in Equation 12 below.

[Equation 12]

는 DBP 데이터 g(t, z)가 있는 입력 공간을 의미하며,

는 영상 f(t, z)가 속한 공간을 의미할 수 있다. Tikhonov 정규화의 경우

는 폐쇄 형태를 가지고 있지만, l₁ 또는 총 변동(TV; total variation)과 같은 일반 정규화 함수에 대해서는 역 매핑

가 일반적으로 비선형이며, 계산적으로 비싼 반복법을 이용하여 찾아야 한다. 런타임 계산 복잡성을 줄이기 위한 빠른 해결책은 비선형 매핑

를 사전 계산하는 것이다. 불행하게도 매핑

는 입력에 의존하기 때문에 모든 입력에 대해

를 저장하는 것은 엄청난 양의 메모리를 필요로 할 것이고 실현 가능하지도 않을 것이다.here,

denotes the input space with DBP data g(t, z),

may mean a space to which the image f(t, z) belongs. For Tikhonov regularization

has a closed form, but is inversely mapped for general regularization functions such as _{l 1 or total variation (TV).}

is generally non-linear and must be found using computationally expensive iterative methods. A quick solution to reduce runtime computational complexity is non-linear mapping

is to be pre-calculated. unfortunately mapping

is input dependent, so for every input

Storing them would require a huge amount of memory and would not be feasible.

를 각 영역에 대한 입력이 공통의 선형 표현을 공유하는 겹치지 않는 영역으로 분할한다. 그런 다음, 각 입력에 대한 해당 선형 표현으로의 전환은 ReLU 활성화 패턴을 기반으로 즉시 수행할 수 있다.In this respect, an encoder-decoder CNN (ED CNN) using rectified linear unit (ReLU) nonlinearity can solve this problem. Specifically, recent theoretical studies (JC Ye and WK Sung, "Understanding geometry of encoder-decoder CNNs," in Proceedings of the 36th International Conference on Machine Learning, ser. Proceedings of In Machine Learning Research, K. Chaudhuri and R. Salakhutdinov, Eds., vol. 97. Long Beach, California, USA: PMLR, 09-15 Jun 2019, pp. 7064-7073.) Applicants of the present invention describe ReLU nonlinear We have shown that encoder-decoder CNNs with gender generate a large number of local linear mappings. More specifically, the input space

split into non-overlapping regions where the inputs to each region share a common linear representation. Then, the conversion to the corresponding linear representation for each input can be performed immediately based on the ReLU activation pattern.

를 추정하는 것으로, 아래 <수학식 13>과 같이 나타낼 수 있다.Another unique aspect of deep neural networks is that these exponentially many linear representations can be derived from a small set of filters thanks to the combinatorial nature of ReLU nonlinearity. In particular, network training is a filter set

, which can be expressed as <Equation 13> below.

[Equation 13]

은 실측 자료(ground truth) 영상과 DBP 영상으로 구성된 트레이닝 데이터 세트를 의미하고,

는

에 의해 매개변수화된 역 매핑을 의미할 수 있다.here,

is a training data set consisting of a ground truth image and a DBP image,

Is

may mean inverse mapping parameterized by .

는 훨씬 더 작은 메모리와 공간을 요구한다. 반면에, 연관된 선형 표현들의 수는 네트워크 깊이, 너비, 스킵된 연결(skipped connection)과 함께 기하급수적으로 증가한다.Filter parameters in general

requires much less memory and space. On the other hand, the number of associated linear representations increases exponentially with network depth, width, and skipped connections.

This representation and input adaptability are considered to be one of the main reasons for the success of deep neural network networks in the image reconstruction problem.

spectral blending (spectral blending)

The factorization method leads to a deep learning approach for each plane of interest, but a previous study (M. Lee, Y. Han, JP Ward, M. Unser, and JC Ye, "Interior tomography using 1d generalized total variation. part ii: Multiscale implementation," SIAM Journal on Imaging Sciences, vol. 8, no. 4, pp. 2452-2486, 2015.) showed that the missing frequency domain depends on the orientation of the plane of interest. However, previous studies also showed that artifacts in the missing frequency domain can be minimized by appropriately combining the missing frequency domains for different filtering directions. This is explained in detail as follows.

Suppose we are interested in recovering voxels from DBP data to x = (x, y, z). 5(a) is a plan view of a point x and a source trajectory. In 4(a), a _v (λ) means a corresponding virtual source position.

가 수평 방향으로 정렬되는 관상 방향(예를 들면, a_v(λ₁ ^-)과 a_v(λ₁ ⁺) 사이)의 영상 재구성을 고려해 본다. 이전 연구(M. Lee, Y. Han, J. P. Ward, M. Unser, and J. C. Ye, "Interior tomography using 1d generalized total variation. part ii: Multiscale implementation," SIAM Journal on Imaging Sciences, vol. 8, no. 4, pp. 2452-2486, 2015.)의 스펙트럴 분석에 따르면, DBP 결과 데이터는 도 5(b)와 같이 시상 방향을 따라 누락된 주파수 영역을 가지고 있다. 따라서, 관상 방향 DBP 데이터를 사용하는 디컨볼루션 알고리즘은 누락된 주파수 영역의 부적절한 정의로 인해 시상 방향을 따라 노이즈 부스팅을 가질 수 있다. 마찬가지로, 시상 방향(예를 들면 a_v(λ₂ ^-)과 a_v(λ₂ ⁺) 사이)의 DBP 데이터에 대해서는 도 5(c)에 해당 주파수 영역이 제시되어 디코볼루션 후 노이즈 부스팅이 발생한다. plane of interest

Consider the image reconstruction in the coronal direction (eg, between a _v (λ ₁ ^- ) and a _v (λ ₁ ^{+ )) in which is horizontally aligned.} A previous study (M. Lee, Y. Han, JP Ward, M. Unser, and JC Ye, "Interior tomography using 1d generalized total variation. part ii: Multiscale implementation," SIAM Journal on Imaging Sciences, vol. 8, no. 4, pp. 2452-2486, 2015.), the DBP result data has a missing frequency domain along the sagittal direction as shown in FIG. Therefore, the deconvolution algorithm using the coronal direction DBP data may have noise boosting along the sagittal direction due to improper definition of the missing frequency domain. Similarly, for DBP data in the sagittal direction (eg, between a _v (λ ₂ ^- ) and a _v (λ ₂ ⁺ )), the corresponding frequency domain is presented in FIG. 5(c), and noise boosting occurs after decovolution. do.

In the case of full-scan cone-beam CT, there is redundancy in the 2D digitization of 3D data. FIG. 5( d ) shows the missing frequency domain by combining the reconstruction results in two directions. In particular, the common missing frequency domain is a dark rectangle centered at the origin, which is the intersection of the two missing frequency domains in coronal and sagittal processing. This means that, except for the common frequency domain, the spectral component missing from coronal direction processing can be compensated for by the sagittal processing result, and vice versa.

6 shows an exemplary diagram for explaining a spectral blending technique, and shows a flowchart of spectral blending to achieve such a synergistic combination. As shown in Fig. 6(a), the axial images (or slice images) f ^cor and f ^sag of the coronal and sagittal reconstructed (or reconstructed) 3D volumes are first transformed into the spectral domain using a 2D Fourier transform. . Since it is improperly defined along the missing frequency domain, the artifact of the straight pattern indicated by the yellow arrow in Fig. 6(b) is usually visible in the Fourier domain along the missing frequency domain. Then, to suppress the signal in the missing frequency domain, a bow-type spectral weight is applied and combined together. This process can be mathematically expressed as <Equation 14> below.

[Equation 14]

및

은 각각 푸리에 변환과 역 푸리에 변환을 의미하고, ω는 보우-타이(bow-tie) 스펙트럴 마스크를 의미하며,

는 요소별 곱(element-wise multiplication)을 의미할 수 있다.here,

and

denotes Fourier transform and inverse Fourier transform, respectively, ω denotes a bow-tie spectral mask,

may mean element-wise multiplication.

data set

In the present invention, a 10-subject data set from the American Association of Phsicists in Medicine (AAPM) low-dose CT Grand Challenge can be used. Numerically obtain cone-beam CT sine wave data using anterior cone-beam projection from volumetric image data. The number of detectors may be ^{1440×1440 elements with a pitch of 1 mm 2} , the number of views may be 1200, the source to source distance (DSO) may be 500 mm, and the source to detector distance (DSD) may be 1000 mm. The maximum cone angle is calculated to be 35.8°. The x - y size of the image is 512Х512 and the z size is 400 to 600, but it may vary from patient to patient. Patient data from 8 out of 10 patients can be used as the training set, 1 patient data can be used as the validation set, and other patient data can be used as the test set. Since the neural network in the present invention is trained in the plane of interest according to coronal slices and sagittal slices, the training data set also consists of coronal slices and sagittal slices. This corresponds to 8192 (=512Х8Х2) and 1024 (=512Х1Х2) coronal/sagittal slices for training and validation, respectively. Of course, the neural network in the present invention is not limited to being trained according to the coronal slice and the sagittal slice, and may be trained according to slices in at least two different directions, and the two directions are not necessarily limited to being orthogonal to each other. In addition, since the neural network in the present invention is trained on slices in two or more directions, it may be trained using slices in three or more directions. Of course, in this case, slices in three or more directions may be used as a neural network input, and slices in two or more directions may be used as a neural network input.

In addition, the present invention may conduct an experiment using a numerical phantom to more quantitatively investigate network performance. 7 shows an exemplary view of a disk phantom, and shows a typical example of a numerical phantom composed of disks having the same thickness and spacing between adjacent disks. The radius of the disk can be 80mm, and various phantoms can be made by changing the size, radius, number, thickness, and spacing as shown in Table 1 below. Here, (a) to (c) of Table 1 may be used for fine nuning, and (d) may be used for testing.

Algorithm implementation

평면의 DBP 영상은 두 개의 보완 소스 궤적, 즉 각 λ[λ^-, λ⁺]와 λ[0, 2π]\[λ^-, λ⁺]에서 a(λ)를 사용하여 얻을 수 있으며, 여기서 \은 차세트, 즉 A\B는 A∩B^c를 나타낸다. 이 경우 도 8과 같이 두 개의 DBP 영상을 모두 입력으로 사용한다.8 shows an exemplary diagram of an input and output of a neural network for removing cone beam artifacts, and shows an exemplary diagram of an input and output of the deep learning method of the present invention. As shown in FIG. 8 , images from two directions, for example, coronal and sagittal view DBP images are used as inputs, and artifact-free coronal and sagittal view images are used as labels. For full scan cone beam CT

A planar DBP image can be obtained using a(λ) from two complementary source trajectories: λ[λ ^- , λ ⁺ ] and λ[0, 2π]\[λ ^- , λ ^{+ ], where \} is the difference set, that is, A\B denotes A∩B ^c . In this case, both DBP images are used as inputs as shown in FIG. 8 .

9 shows an exemplary diagram of a deep neural network architecture in the method of the present invention, and the same neural network, for example, one neural network may be used for each DBP image, that is, the initial image.

The network backbone corresponds to the modified architecture of U-Net. The neural network in FIG. 9 is a convolution layer that performs a linear transform operation, a batch normalization layer that performs a normalization operation, and a ReLU that performs a nonlinear function operation. Includes a (rectified linear unit) layer and a contracting path connection with concatenation. In particular, each stage contains 4 sequential layers consisting of Convolution, Batch Normalization and ReLU layers with 3 Х 3 kernels. The last stage includes two sequential layers and a last layer, and the last layer includes only a convolution with a 1 Х 1 kernel that generates a reconstructed image for the first image. The number of channels for each convolutional layer is shown in Fig. 9, the number of channels is doubled after each pooling layer, and the size of the layers is reduced by a factor of four. Here, the pooling layer may be a 2 Х 2 average pooling layer, and the pooling layer may be located between each stage. Also, the unpooling layer may be a 2 Х 2 average unpooling layer. The total number of trainable parameters is about 22,000,000.

The present invention can be compared by implementing two additional algorithms, for example, a deep convolutional neural network using FDK reconstruction (FDK-CNN) and an MBIR algorithm to which a TV penalty is imposed. FDK-CNN was trained to learn artifact-free images from the reconstructed images of the FDK algorithm using incomplete projection data due to cone angle. The input image is an FDK image damaged by cone-beam artifacts, while the artifact-free data is used as the actual measurement data. For fair comparison, coronal and sagittal view images are used as network inputs similar to the DBP domain network of the present invention. Since FDK reconstruction combines projections from all full scan angles, it is not possible to obtain two different FDK images for each processing direction. This is why only two images along the coronal and sagittal views are provided as inputs. The same U-Net architecture can be used for FDK domain neural networks. Therefore, the only difference from the DBP domain deep network of the present invention is the input image.

Reconstruction to which a total variation (TV) penalty is imposed can be implemented by solving an optimization problem as shown in Equation 15 below.

[Equation 15]

Here, y and f mean the measured sinogram data and its 3D image, A means the system matrix, ∇ _{x, y, z} means the difference operator along the (x, y, z) axis, , λ may mean a normalization factor.

For quantitative evaluation, a normalized mean square error (NMSE) value may be used, and this value may be expressed as in Equation 16 below.

[Equation 16]

Here, f and f ^* may mean a reconstructed image and a measured image, respectively, and M and N may mean the number of rows and columns of pixels.

Also, a peak signal-to-noise ratio (PSNR) defined as in Equation 17 below may be used.

[Equation 17]

In addition, a structural similarity index defined as in Equation 18 below may be used.

[Equation 18]

where μ _f means the mean of f, σ ² _f means the variance of f, and σ _ff* is means the cross covariance of f and f*, the variables for stabilizing the distribution are c ₁ = (k ₁ L) ² and c ₂ = (k ₂ L) ² , L means the dynamic range of pixel intensity, k ₁ and k ₂ are constants, and the default values may be k ₁ = 0.01 and k ₂ = 0.03.

network training

The present invention can use MatConvNet toolbox (ver.24) to implement FDK and DBP domain networks in MATLAB R2015a environment, and the processing unit is an Intel Core i7-7700 (3.60GHz) central processing unit (CPU) and GTX 1080Ti graphics It may be a processing unit (GPU). l ₂ The loss may be used as an objective function of training, and the number of training network epochs may be 300. Stochastic gradient descent (SGD) can be used as an optimizer for network training. The initial learning rate may be 10 ^{-4 , and may gradually drop to 10 -5} at each epoch, and the regularization parameter may be ^{10 -4.} For data augmentation, the entire data set can be performed with horizontal and vertical flips. The mini-batch is used as 4, and the size of the input patch is 256Х256. Since the convolution operation is spatially invariant, the trained filter can be used on the entire input data in the estimation step. The training time may be a period of time, for example, about 4 days.

10 shows an exemplary diagram of the convergence of the objective function according to the number of epochs, wherein the dotted line and the solid line indicate the objective function of the training and verification steps, respectively. As shown in Figure 10, it can be seen that the network of the present invention is well trained and not over-fitted because the training and validation curves converge closely.

Since numerical phantoms are different from AAPM data sets, FDK-CNN and our network can be fine-tuned using subsets of disk phantoms to evaluate the reconstruction performance of other disk phantoms. Specifically, the disk phantom using the parameters of (a) (b) (c) of Table 1 can be used for precise adjustment, and the disk phantom of (d) can be used in the test step.

11 and 12 show exemplary views comparing restoration performance by iterative restoration technique (total variation), FDK-CNN and the method of the present invention, and various quantitative indicators such as restoration performance and PSNR, NMSE and SSIM. will be.

11 and 12, in the case of the FDK reconstruction image, it can be seen that the image is darkened and the shape is distorted in the portion having a large cone angle. In the case of the TV restoration method, the intensity is restored, but the It can be seen that the texture and morphological distortions are not sufficiently restored. In the case of FDK-CNN, the intensity and texture are restored, but the morphological distortion is not restored, while the method of the present invention completely restores the intensity and texture as well as the distortion of the image.

In addition, as can be seen from the quantitative indicators in <Table 2> below, it can be seen that the method of the present invention shows excellent performance in all quantitative indicators compared to other algorithms.

In addition, the present invention can use a neural network based on a convolutional framelet, and the convolutional framelet is a method of expressing an input signal using a local basis and a non-local basis, which will be described as follows.

)와 비국소 기저 (

)를 이용하여 표현한 것으로, 아래 <수학식 19>와 같이 나타낼 수 있다.The convolutional framelet is a local basis (

) and non-localized basal (

), which can be expressed as in <Equation 19> below.

[Equation 19]

는 비국소 기저 벡터를 가지는 선형 변환 연산을 의미하고,

는 국소 기저 벡터를 가지는 선형 변환 연산을 의미할 수 있다.here,

is a linear transformation operation with a non-local basis vector,

may mean a linear transformation operation having a local basis vector.

와

를 가질 수 있으며, 기저 벡터들의 직교 관계는 아래 <수학식 20>과 같이 정의될 수 있다.In this case, the local basis vector and the non-local basis vector are dual basis vectors orthogonal to each other, respectively.

Wow

, and the orthogonal relationship between the basis vectors may be defined as in Equation 20 below.

[Equation 20]

Using Equation 20, the convolutional framelet can be expressed as in Equation 21 below.

[Equation 21]

는 행켈 행렬 연산(Hankel matrix operator)을 의미하는 것으로, 컨볼루션 연산을 행렬곱(matrix multiplication)으로 표현할 수 있게 해주며, C는 국소 기저와 비국소 기저에 의하여 변환된 신호인 컨볼루션 프레임렛 계수(convolution framelet coefficient)를 의미할 수 있다.here,

stands for the Hankel matrix operator, which allows the convolution operation to be expressed as matrix multiplication, and C is the convolutional framelet coefficient, which is a signal transformed by a local basis and a non-local basis. (convolution framelet coefficient).

를 적용하여 본래의 신호로 복원될 수 있다. 신호 복원 과정은 아래 <수학식 22>와 같이 나타낼 수 있다.Convolutional framelet coefficients C are dual basis vectors

can be restored to the original signal by applying . The signal restoration process can be expressed as in Equation 22 below.

[Equation 22]

As described above, a method of representing an input signal through a local basis and a non-local basis is called a convolutional framelet.

As described above, the method according to an embodiment of the present invention removes the cone beam artifact with a cone beam tomography image using a neural network that removes the cone beam artifact, and uses a spectral blending technique in which the cone beam image from which the cone beam artifact is removed is input. By removing the residual noise, the cone-beam tomography image can be restored to a high-quality image.

13 is a diagram illustrating a configuration of an image processing apparatus according to an embodiment of the present invention, and shows a configuration of an apparatus performing the method of FIGS. 1 to 12 .

Referring to FIG. 13 , an image processing apparatus 1300 according to an embodiment of the present invention includes a receiving unit 1310 and a restoration unit 1320 .

The receiver 1310 receives cone beam computed tomography data, for example, projection data.

The restoration unit 1320 reconstructs an image of the cone-beam computed tomography data using a neural network that removes cone-beam artifacts.

At this time, the restoration unit 1320 sets the internal computed tomography data to the initial images in at least two directions from the cone beam computed tomography data using a preset analytic algorithm, for example, the coronal direction and the sagittal direction. First images for each (sagittal) direction are acquired, and an image (CNN output) from which cone beam artifacts are removed from the first images for at least two directions using a neural network inputting the first images (CNN input) By outputting the image, it is possible to reconstruct the image for the cone-beam computed tomography data. Here, the restoration unit 1320 may acquire initial images in two orthogonal directions or may acquire initial images in two non-orthogonal directions, and by using a neural network to remove cone beam artifacts in each of the initial images. , reconstructed images in each of two directions from which cone beam artifacts are removed may be obtained. The analytic algorithm may include a differential backprojection (DBP) algorithm.

The neural network used by the reconstruction unit 1320 may include a convolution framelet-based neural network, and may include a multi-resolution neural network including a pooling layer and an unpooling layer. can In addition, the neural network may include a bypass connection from the pooling layer to the unpooling layer.

Furthermore, the reconstruction unit 1320 generates a final reconstructed image for the reconstructed image using the spectral blending technique.

Here, the restoration unit 1320 reconstructs a 3D image through initial images for each direction by using each of the initial images for at least two directions from which the cone beam artifact has been removed, and another 3D image for each of the reconstructed 3D images. A final reconstructed image may be generated by removing residual noise using a spectral blending technique that takes a slice image pair in the direction as an input.

Although the description of the device of FIG. 13 is omitted, each component constituting FIG. 10 may include all the contents described with reference to FIGS. 1 to 12 , which is apparent to those skilled in the art.

The device described above may be implemented as a hardware component, a software component, and/or a combination of the hardware component and the software component. For example, devices and components described in the embodiments may include, for example, a processor, a controller, an arithmetic logic unit (ALU), a digital signal processor, a microcomputer, a field programmable array (FPA), It may be implemented using one or more general purpose or special purpose computers, such as a programmable logic unit (PLU), microprocessor, or any other device capable of executing and responding to instructions. The processing device may execute an operating system (OS) and one or more software applications running on the operating system. The processing device may also access, store, manipulate, process, and generate data in response to execution of the software. For convenience of understanding, although one processing device is sometimes described as being used, one of ordinary skill in the art will recognize that the processing device includes a plurality of processing elements and/or a plurality of types of processing elements. It can be seen that can include For example, the processing device may include a plurality of processors or one processor and one controller. Other processing configurations are also possible, such as parallel processors.

The software may comprise a computer program, code, instructions, or a combination of one or more thereof, which configures a processing device to operate as desired or is independently or collectively processed You can command the device. The software and/or data may be any kind of machine, component, physical device, virtual equipment, computer storage medium or device, to be interpreted by or to provide instructions or data to the processing device. may be embodied in The software may be distributed over networked computer systems, and stored or executed in a distributed manner. Software and data may be stored in one or more computer-readable recording media.

The method according to the embodiment may be implemented in the form of program instructions that can be executed through various computer means and recorded in a computer-readable medium. The computer-readable medium may include program instructions, data files, data structures, etc. alone or in combination. The program instructions recorded on the medium may be specially designed and configured for the embodiment, or may be known and available to those skilled in the art of computer software. Examples of the computer-readable recording medium include magnetic media such as hard disks, floppy disks and magnetic tapes, optical media such as CD-ROMs and DVDs, and magnetic such as floppy disks. - includes magneto-optical media, and hardware devices specially configured to store and carry out program instructions, such as ROM, RAM, flash memory, and the like. Examples of program instructions include not only machine language codes such as those generated by a compiler, but also high-level language codes that can be executed by a computer using an interpreter or the like.

As described above, although the embodiments have been described with reference to the limited embodiments and drawings, various modifications and variations are possible by those skilled in the art from the above description. For example, the described techniques are performed in a different order than the described method, and/or the described components of the system, structure, apparatus, circuit, etc., are combined or combined in a different form than the described method, or other components Or substituted or substituted by equivalents may achieve an appropriate result.

Therefore, other implementations, other embodiments, and equivalents to the claims are also within the scope of the following claims.

Claims (16)

Receiving the cone beam computed tomography data in the receiver; and
reconstructing an image of the cone beam computed tomography data using a neural network that removes cone beam artifacts in a restoration unit
includes,
Restoring the image
An initial image in a first direction and an initial image in a second direction orthogonal to the first direction are obtained using a differential inverse projection algorithm, and the neural network is performed with respect to the initial image in the first direction and the initial image in the second direction. reconstructing a first image from which the cone beam artifact in the first direction is removed and a second image from which the cone beam artifact in the second direction is removed using a network;
Restoring the image
A first image in which the cone beam artifact in the first direction is removed using a spectral blending technique in which a first image from which the cone beam artifact in the first direction is removed and a second image from which the cone beam artifact in the second direction is removed are input Removes noise remaining in the first image and the second image from which the cone beam artifacts in the second direction are removed, and slices in the third direction using the image in the first direction and the image in the second direction from which the residual noise is removed An image processing method of acquiring images and generating a final reconstructed image of the cone beam computed tomography data using the obtained slice images in the third direction.
delete delete delete delete According to claim 1,
The neural network is
An image processing method comprising a neural network based on a convolution framelet.
According to claim 1,
The neural network is
An image processing method comprising a multi-resolution neural network including a pooling layer and an unpooling layer.
8. The method of claim 7,
The neural network is
and bypass connection from the pooling layer to the unpooling layer.
a receiver for receiving cone-beam computed tomography data; and
A restoration unit that restores an image of the cone-beam computed tomography data using a neural network that removes cone-beam artifacts
includes,
the restoration unit
An initial image in a first direction and an initial image in a second direction orthogonal to the first direction are obtained using a differential inverse projection algorithm, and the neural network is performed with respect to the initial image in the first direction and the initial image in the second direction. reconstructing a first image from which the cone beam artifact in the first direction is removed and a second image from which the cone beam artifact in the second direction is removed using a network;
the restoration unit
A first image in which the cone beam artifact in the first direction is removed using a spectral blending technique in which a first image from which the cone beam artifact in the first direction is removed and a second image from which the cone beam artifact in the second direction is removed are input Removes noise remaining in the first image and the second image from which the cone beam artifacts in the second direction are removed, and slices in the third direction using the image in the first direction and the image in the second direction from which the residual noise is removed An image processing apparatus for acquiring images and generating a final reconstructed image of the cone-beam computed tomography data using the obtained slice images in the third direction.
delete delete delete delete 10. The method of claim 9,
The neural network is
An image processing apparatus comprising a neural network based on a convolution framelet.
10. The method of claim 9,
The neural network is
An image processing apparatus comprising a multi-resolution neural network including a pooling layer and an unpooling layer.
16. The method of claim 15,
The neural network is
and bypass connection from the pooling layer to the unpooling layer.

Images