CN109035156B - DNST (deep depth transform) -based medical CT (computed tomography) image denoising method - Google Patents

DNST (deep depth transform) -based medical CT (computed tomography) image denoising method Download PDF

Info

Publication number
CN109035156B
CN109035156B CN201810618315.0A CN201810618315A CN109035156B CN 109035156 B CN109035156 B CN 109035156B CN 201810618315 A CN201810618315 A CN 201810618315A CN 109035156 B CN109035156 B CN 109035156B
Authority
CN
China
Prior art keywords
image
shear wave
coefficient
noise
dnst
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810618315.0A
Other languages
Chinese (zh)
Other versions
CN109035156A (en
Inventor
程芸
周俊
田峥
赵恺伦
张聚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang Hospital
Original Assignee
Zhejiang Hospital
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang Hospital filed Critical Zhejiang Hospital
Priority to CN201810618315.0A priority Critical patent/CN109035156B/en
Publication of CN109035156A publication Critical patent/CN109035156A/en
Application granted granted Critical
Publication of CN109035156B publication Critical patent/CN109035156B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20036—Morphological image processing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

The DNST-based medical CT image denoising method comprises the following steps: step 1) establishing a medical CT image model; step 2) decomposing the CT image by dispersing the inseparable shear wave; step 3) carrying out bivariate contraction treatment on the high-frequency shear wave coefficient; step 4) carrying out rotation invariant bilateral non-local mean filtering processing on the low-frequency shear wave coefficient, and step 4) carrying out inverse DNST (discrete Fourier transform) on the processed coefficient. The invention is compared with the latest denoising field algorithm through experimental analysis, and is effectively applied to the medical CT denoising field; particularly, rotation invariant bilateral non-local mean filtering and bivariate contraction of high-frequency shear wave coefficients are adopted for CT low frequency, and detail information of soft tissues in medical CT images can be better protected. Through comparison of a large amount of experimental data, a medical CT image denoising method based on combination of shear wave domain bivariate shrinkage and rotation invariant bilateral non-local mean filtering is provided, and analysis and diagnosis of doctors can be better facilitated.

Description

DNST (deep depth transform) -based medical CT (computed tomography) image denoising method
Technical Field
The invention relates to the field of medical image denoising, in particular to a medical CT image. A medical CT image denoising method based on shear wave domain bivariate shrinkage and rotation invariant bilateral non-local mean filtering combination is designed.
Background
In the history of human development, medical science has been the subject of much attention and continuous development. The popularity and development of medical image digitization have greatly improved the efficiency and accuracy of medical diagnosis. Among them, Computed Tomography (CT) has a strong capability in radiation diagnosis, and it is the first one to obtain an axial slice of the internal structure of a human body without opening the human body, and it has the advantages of high density resolution of tissue structure, little damage to the human body, and easy operation, and is widely used in the diagnosis of cancer, cardiovascular disease, infectious disease, appendicitis, trauma and musculoskeletal diseases, etc. The medical CT discovery created the importance of three-dimensional nobel prize reflecting CT imaging technology from the side.
CT is a medical imaging method using tomography, which can simply acquire internal structural features of an object, such as size, shape, internal defects and density, and is a powerful non-destructive evaluation technique widely used to generate two-dimensional and three-dimensional cross-sectional images of the human body. CT imaging uses physical technology, and fully utilizes the characteristic that X-rays are absorbed in a human body on the basis of measuring the attenuation coefficient of the X-rays in the human body. The received ray is converted into electric signal by detector in medical instrument, and the electric signal is converted into ray attenuation coefficient value by analog-digital conversion. And (3) solving the gray value of an image picture of a certain section of the human body by a computer according to the obtained two-dimensional distribution matrix of the ray attenuation coefficient value by adopting an excellent mathematical algorithm so as to obtain a sectional view of a certain part of the human body. However, the accuracy of CT positioning and qualitative diagnosis still receives the influence of the lesion site, size, nature, length of the course of disease, and the like, and thus the functional information of the visceral organs cannot be reflected. The image is also easily affected by factors such as noise interference, shot details and the like, so that spots and fine grains are observed in the image, and further, the doctor is influenced to distinguish certain size details of the focus part from the background. Noise interference in the image will affect further analysis of the image, so that the noise is filtered from the CT image by researching a related medical image denoising algorithm, and a good visual environment can be provided for observing the image.
The present invention uses medical CT as the subject of study. Since the CT imaging is not immune to various physical factors, the quality of the CT image is seriously affected by the existence of speckle noise, resulting in poor medical image quality. The presence of noise masks image features that have little difference in gray scale. Noise can be a significant distraction to clinicians for their accurate diagnosis, especially for less experienced physicians. Therefore, from the perspective of clinical application, it is necessary to research a CT medical image denoising method, so as to provide technical support for doctors to make more accurate diagnosis and reduce the risk of manual diagnosis.
In conclusion, the method for denoising the medical CT image has wide application.
Disclosure of Invention
The invention provides a denoising algorithm for a discrete inseparable shear wave (DNST) medical CT image, which aims to overcome the defects of the sparse capability of wavelet analysis processing high-dimensional data and the poor boundary effect of a discrete shear wave frame and is used for denoising the medical CT image.
In the prior art, wavelet change can be well used for image denoising and effectively seizes one-dimensional singular points, but cannot reflect abrupt changes of straight lines and curves. The ridge wave transformation can well capture the singularity of the line, make up for the deficiency of the wavelet, but still can not effectively capture the singularity of the curve. In recent years, medical images are processed through a discrete shear wave algorithm, so that the medical image denoising technology field is broken through to a certain extent. DNST has better frame bounds and better directional selection than discrete separable shear waves, meaning DNST has better denoising effect. The DNST tool kit is applied to medical CT noise images, the traditional shrinkage algorithm is modified to use a bivariate shrinkage algorithm, and rotation-invariant bilateral non-local mean filtering is performed in the low-frequency part combined with the CT images. The invention designs a medical CT image denoising method based on the combination of shear wave domain bivariate shrinkage and rotation invariant bilateral non-local mean filtering, which has the advantages of high speed and obvious denoising, and finally verifies the feasibility and the optimization effect of the method through simulation.
Compared with the prior art, the invention has the creativity and advantages that: a medical CT image denoising method based on shear wave domain bivariate shrinkage and rotation invariant bilateral non-local mean filtering combination is provided, and DNST transformation overcomes the defects of the sparse capability of wavelet analysis processing high-dimensional data and the boundary of a discrete shear wave frame. Compared with the traditional hard threshold and soft threshold, the double variable shrinkage algorithm is more suitable for being applied to the field of CT denoising, and rotation invariant bilateral non-local mean filtering suitable for CT characteristics is added to the low-frequency part. The method has the advantages of multi-resolution, multi-scale, multi-directivity and time-frequency locality, can better protect the image edge information when being applied to the CT image denoising, and provides convenience for the diagnosis of doctors.
In order to make the objects, technical solutions and advantages of the present invention more clear, the technical solutions of the present invention are further described below and are divided into the following four steps.
The DNST-based medical CT image denoising method comprises the following steps:
step 1) establishing a medical CT image model;
in order to solve the problem of CT noise, the user cannot judge the CT noise by means of subjective feeling of the user, and the CT noise can be recognized only by a probability statistical method. Therefore, the most important is to embody the abstract noise in the CT image and to build a mathematical model conforming to the basic characteristics.
The energy of the X-rays is transmitted in the form of a single energy block of quanta, so that these infinite X-ray quanta are detected by the X-ray detector. Due to statistical errors, the detected X-ray quanta may differ for one more measurement. Statistical noise in CT images arises from fluctuations that may occur in the detection of small numbers of X-ray quantum dots, and thus may also be referred to as quantum noise. This is gaussian additive noise in terms of a mathematical model, which is as follows:
Inoisy(x,y)＝uorignal(x,y)+ηnoisy(x,y) (1)
where (x, y) denotes a pixel point in the image, InosiyRepresenting CT images containing noise, i.e. observed CT images, uorignalRepresenting noiseless CT images, etanoisyRepresenting noise, the distribution signAnd a Gaussian distribution.
Step 2) decomposing the CT image by dispersing the inseparable shear wave;
s1: inputting a two-dimensional CT image signal f ∈ RX*YA scale parameter J ∈ N, a direction parameter k ∈ NJAnd selecting a directional filter DirectionFilter and a low-pass filter quadraturecrrorfilter.
S2: converting the frequency spectrum f of a CT imagefreq＝FFT(f)。
S3: calculating sub-band i epsilon [0, nth ] under different scales and different directions]Shear wave coefficient under (i) sheearlettCoeffs (i) epsilon RX*Y*nthAccording to the convolution theory of the shear wave domain, the following can be obtained:
Figure BDA0001697446790000031
s4: the discrete unseparatable shear wave coefficient shearlet coeffs (i) is obtained by calculation.
Wherein nth in step 3 represents the redundancy of the entire tightly supported DNST system, which is calculated as follows:
nth＝2*((2*2k[0]+1))+2*((2*2k[1]+1))+...+2*((2*2k[J]+1) (3)
DNST enables better frame definition. In addition to this, shear waves generated by an inseparable shear wave generator
Figure BDA0001697446790000032
One major advantage of (a) is that the directional selectivity is improved for each size of the sector filter P in the frequency domain. Decomposition of medical CT images into f in the frequency domain by DNST decomposition1,f2,...,fnth-1High-frequency CT image and low-frequency CT image f with equal sheet sizenth。
Step 3) on the high frequency shear wave coefficient f1,f2,...,fnth-1Carrying out bivariate shrinkage treatment;
one of the strongest techniques for image denoising is a bivariate threshold function denoising method commonly used in wavelet domain. Here, the shear wave domain is extended with this idea. The shear wave coefficient model is calculated by using a non-gaussian bivariate function. The relationship of the finer and coarse scales is first defined and then the shear wave coefficients are modeled efficiently, each coefficient depending on a coefficient at the same spatial position in the coarser scale, called CS, and each CS depending on a coefficient in the same spatial position in the immediate finer scale, called FS. Thus, there is one CS per CS in the shear domain and four FS per CS. The model can effectively capture the dependency between the frame coefficients and their CS.
In general, for a given original noise-free image, if it is disturbed by additive white gaussian noise, the degraded image is represented in the shear domain as follows:
g＝f+ε (4)
where g, f, epsilon represent the observed shear wave coefficient, the original noise-free shear wave coefficient, and the noise coefficient, respectively. The purpose of denoising is to obtain an estimate of the original coefficients
Figure BDA0001697446790000041
So that the estimated value
Figure BDA0001697446790000042
As close as possible to the original noise-free coefficient f. In the proposed method, maximum a posteriori probability (MAP) is used to estimate the de-noising coefficients
Figure BDA0001697446790000043
Therefore:
Figure BDA0001697446790000044
bayesian rules can be used to write:
Figure BDA0001697446790000045
the following equation can thus be derived:
Figure BDA0001697446790000046
to describe the dependency between CS and FS in the framework domain, let f1Denotes f2CS (note f)1Is and f2Frame coefficients at the same spatial location, but at the next finer scale), so there is f ═ f (f ═ f)1,f2)，g＝(g1,g2)，ε＝(ε1,ε2)，f1And f2Are not relevant, but are not independent. Thus g ═ f + epsilon can be written in the form:
Figure BDA0001697446790000047
assume that the mean and variance of the additive noise are 0 and
Figure BDA0001697446790000048
from pεThe noise probability density function represented by (ε) is as follows:
Figure BDA0001697446790000049
from the distribution of the shear wave coefficients, it is possible to pass through the bivariate probability density function p from a symmetrical circular probability distributionf(f) And (6) fitting the model. Thus, it can be expressed as follows:
Figure BDA00016974467900000410
substituting equations (9) and (10) into equation (7) yields:
Figure BDA00016974467900000411
if p isf(f) Is convex, solves forThe process equation (11) can be converted to the following solving equation:
Figure BDA0001697446790000051
wherein
Figure BDA0001697446790000052
Substituting equation (13) into equation (12), f1The MAP estimator of (a) is a bivariate contraction function:
Figure BDA0001697446790000053
from this equation, it can be seen that for each shear coefficient there is a corresponding threshold value, which depends not only on the CS coefficient but also on the FS coefficient.
Step 4) on the low frequency shear wave coefficient fnthCarrying out rotation invariant bilateral non-local mean filtering processing;
in 2009, Coipe et al proposed an Optimized Bayesian non-local Means (OBNLM) filter for speckle noise reduction of images based on block-by-block implementation of CNLM using Pearson distance metric. However, OBNLM filters are not effective in removing noise from images containing fine structural details and result in excessive smoothing of edges and textures present in the image. To overcome the above limitations, we propose to use a Rotation Invariant Bilateral non-local mean Filter (RIBNLM) proposed by Manjon et al as a post-processing step that optimizes the trade-off between the required smoothness of the individual regions and preserves the details of the image. The proposed filter uses a rotation invariant similarity measure based on intensity values and corresponding neighborhoods to supplement the mean value as
Figure BDA00016974467900000510
And
Figure BDA00016974467900000511
as follows:
Figure BDA0001697446790000054
Figure BDA0001697446790000055
wherein SW represents a search window of 11 × 11 pixels in size, wRIBNLM(i, j) is a weight, and the condition 0. ltoreq. w is satisfiedRIBNLM(i, j) is less than or equal to 1 and sigma-wRIBNLM(i,j)＝1，
Figure BDA0001697446790000056
And
Figure BDA0001697446790000057
is the mean of the reference blob and the blob being processed in the search window SW, and the parameter r is the smoothing parameter. Local mean
Figure BDA0001697446790000058
And
Figure BDA0001697446790000059
the entire image is computed once with only 5 x 5 neighborhoods and stored in an array with the size of the image. Computing using a small neighborhood
Figure BDA0001697446790000061
And
Figure BDA0001697446790000062
to prevent excessive smoothing of the image's singular structures. The RIBNLM filter has a low computational complexity and gives appropriate weights to the blobs, which are structurally similar but have different orientations relative to the reference blob.
Step 5) carrying out inverse DNST conversion on the processed coefficients;
the high-frequency sub-band and the low-frequency sub-band after threshold shrinkage in the above steps are subjected to inverse DNST transformation, so that a denoised CT image which is beneficial to analysis of a doctor can be obtained, and the superiority of the algorithm of the invention in denoising the medical CT image is verified by comparing experimental data. The specific algorithm process of the DNST inverse process is described as follows:
s1 inputting the shear coefficient shearlet coeffs (i) epsilon R after DNST processingX*Y*nth；
S2 setting frec∈RX*YRepresenting the reconstructed image sequence;
s3, calculating each index i to be 0, nth]Reconstructed image sequence frequency spectrum f of lower shearlet coeffs (i)recAnd sum frecAccording to the convolution theory and the frame theory
Figure BDA0001697446790000063
S4, inverse transform of DNST is carried out to obtain a reconstructed image sequence frec:＝IFFT(frec)；
The invention has the following advantages:
1. the invention uses the discrete inseparable shear wave transformation model, has better frame boundary and direction selectivity, and can better effectively and properly decompose the medical CT image.
2. The invention carries out rotation invariant bilateral non-local mean filtering processing on the specific white Gaussian noise in the low-frequency part of the medical CT.
3. In the invention, a self-adaptive threshold algorithm and a bivariate shrinkage algorithm are adopted for the high-frequency part, the relation between coefficients is fully utilized, and a proper shear wave coefficient can be well removed.
4. The method has definite steps and simple structure and has perfect theoretical support.
Drawings
FIG. 1 is a CT in which a noisy image is the sum of a noiseless CT image and a noise component;
FIG. 2 is a flow chart of discrete non-separable shear decomposition;
FIG. 3 is a bivariate shrinkage function model;
FIG. 4 is a schematic diagram of a rotation invariant bilateral non-local mean filtering;
FIG. 5 is a flowchart illustrating the overall steps of the present invention;
fig. 6a to 6f are comparison of experimental results of various algorithms on the CT (σ ═ 40) head, where fig. 6a is a CT noise map, fig. 6b is an NSST algorithm effect map, fig. 6c is an SWT algorithm effect map, fig. 6d is an FDCT algorithm effect map, fig. 6e is an FFST algorithm effect map, and fig. 6f is an algorithm effect map of the present invention;
fig. 7 is a comparison of the denoising effect of the local region.
Detailed Description
The invention is further described below with reference to the accompanying drawings.
The invention discloses a DNST-based medical CT image denoising method, which comprises the following steps:
step 1) establishing a medical CT image model
In order to solve the problem of CT noise, the user cannot judge the CT noise by means of subjective feeling of the user, and the CT noise can be recognized only by a probability statistical method. Therefore, the most important is to embody the abstract noise in the CT image and to build a mathematical model conforming to the basic characteristics.
The energy of the X-rays is transmitted in the form of a single energy block of quanta, so that these infinite X-ray quanta are detected by the X-ray detector. Due to statistical errors, the detected X-ray quanta may differ for one more measurement. Statistical noise in CT images arises from fluctuations that may occur in the detection of small numbers of X-ray quantum dots, and thus may also be referred to as quantum noise. This is gaussian additive noise in terms of a mathematical model, which is as follows:
Inoisy(x,y)＝uorignal(x,y)+ηnoisy(x,y) (1)
where (x, y) denotes a pixel point in the image, InosiyRepresenting CT images containing noise, uorignalRepresenting noiseless CT images, etanoisyRepresenting noise, the distribution of which corresponds to GaussianAnd (3) cloth. The model is shown in fig. 1, and the observed CT image is the sum of a clean CT image and gaussian-distributed noise.
Step 2) decomposition of CT images by means of discrete inseparable shear waves
S1: inputting a two-dimensional CT image signal f ∈ RX*YThe scale parameter J belongs to N, and a shearing vector parameter k belongs to NJAnd selecting a directional filter DirectionFilter and a low-pass filter quadraturecrrorfilter.
S2: calculating a frequency spectrum f of an input signalfreq＝FFT(f)。
S3: calculating sub-band i e [0, nth ] of different scales]Shear wave forward transform coefficients shearlet coeffs (i) e RX*Y*nthAccording to the convolution theory and the framework theory:
Figure BDA0001697446790000071
s4: the discrete unseparatable shear wave coefficient shearlet coeffs (i) is obtained by calculation.
Wherein nth in step 3 represents the redundancy of the entire tightly supported DNST system, which is calculated as follows:
nth＝2*((2*2k[0]+1))+2*((2*2k[1]+1))+...+2*((2*2k[J]+1) (3)
DNST enables better frame definition. In addition to this, shear waves generated by an inseparable shear wave generator
Figure BDA0001697446790000087
One major advantage of (a) is that the directional selectivity is improved for each size of the sector filter P in the frequency domain. As shown in FIG. 2, the medical CT image can be decomposed into f in the frequency domain by continuously decomposing the low frequency part through the discrete inseparable shear wave decomposition1,f2,...,fnth-1High-frequency CT image and low-frequency CT image f with equal sheet sizenth。
Step 3) on the high frequency shear wave coefficient f1,f2,...,fnth-1Carry out bivariateShrink treatment
One of the strongest techniques for image denoising is a bivariate threshold function denoising method commonly used in wavelet domain. Here, the shear wave domain is extended with this idea. The shear wave coefficient model is calculated by using a non-gaussian bivariate function. The relationship of the finer and coarse scales is first defined and then the shear wave coefficients are modeled efficiently, each coefficient depending on a coefficient at the same spatial position in the instantaneous coarser scale, called CS, and each CS depending on a coefficient in the same spatial position in the instantaneous finer scale, called FS. Thus, there is one CS per CS in the shear domain and four FS per CS. The model can effectively capture the dependency between the frame coefficients and their CS.
In general, for a given original noise-free image, if it is disturbed by additive white gaussian noise, the degraded image is represented in the shear domain as follows:
g＝f+ε (4)
where g, f, epsilon represent the observed shear wave coefficient, the original noise-free shear wave coefficient, and the noise coefficient, respectively. The purpose of denoising is to obtain an estimate of the original coefficients
Figure BDA0001697446790000081
So that the estimated value
Figure BDA0001697446790000082
As close as possible to the original noise-free coefficient f. In the proposed method, maximum a posteriori probability (MAP) is used to estimate the de-noising coefficients
Figure BDA0001697446790000083
Therefore:
Figure BDA0001697446790000084
bayesian rules can be used to write:
Figure BDA0001697446790000085
the following equation can thus be derived:
Figure BDA0001697446790000086
to describe the dependency between CS and FS in the framework domain, let f1Denotes f2CS (note f)1Is and f2Frame coefficients at the same spatial location, but at the next finer scale), so there is f ═ f (f ═ f)1,f2)，g＝(g1,g2)，ε＝(ε1,ε2)，f1And f2Are not relevant, but are not independent. Thus g ═ f + epsilon can be written in the form:
Figure BDA0001697446790000091
assume that the mean and variance of the additive noise are 0 and
Figure BDA0001697446790000092
from pεThe noise probability density function represented by (ε) is as follows:
Figure BDA0001697446790000093
from the distribution of the shear wave coefficients, it is possible to pass through the bivariate probability density function p from a symmetrical circular probability distributionf(f) And (6) fitting the model. Thus, it can be expressed as follows:
Figure BDA0001697446790000094
substituting equations (9) and (10) into equation (7) yields:
Figure BDA0001697446790000095
if p isf(f) Being convex, the solution process equation (11) can be converted to the following solution equation:
Figure BDA0001697446790000096
wherein
Figure BDA0001697446790000097
Substituting equation (13) into equation (12), f1The MAP estimator of (a) is a bivariate contraction function:
Figure BDA0001697446790000098
from this equation, it can be seen that for each shear coefficient there is a corresponding threshold, as shown in the bivariate contracting function model of fig. 3, which depends not only on the CS coefficient but also on the FS coefficient.
Step 4) on the low frequency shear wave coefficient fnthPerforming rotation invariant bilateral non-local mean filtering
In 2009, Coipe et al proposed an optimized Bayesian non-local means (OBNLM) filter for speckle noise reduction of images based on block-by-block implementation of CNLM using Pearson distance metric. However, OBNLM filters are not effective in removing noise from images containing fine structural details and result in excessive smoothing of edges and textures present in the image. To overcome the above limitations, we propose to use a Rotation Invariant Bilateral non-local mean Filter (RIBNLM) proposed by Manjon et al as a post-processing step that optimizes the trade-off between the required smoothness of the individual regions and preserves the details of the image. The proposed filter uses intensity-based values and correspondingThe effective rotation-invariant similarity measure of the neighborhood is complemented by the mean value
Figure BDA0001697446790000101
And
Figure BDA0001697446790000102
as follows:
Figure BDA0001697446790000103
Figure BDA0001697446790000104
wherein SW represents a search window of 11 × 11 pixels in size, wRIBNLM(i, j) is a weight, and the condition 0. ltoreq. w is satisfiedRIBNLM(i, j) is less than or equal to 1 and sigma-wRIBNLM(i,j)＝1，
Figure BDA0001697446790000105
And
Figure BDA0001697446790000106
is the mean of the reference blob and the blob being processed in the search window SW, and the parameter r is the smoothing parameter. Local mean as shown in FIG. 4
Figure BDA0001697446790000107
And
Figure BDA0001697446790000108
the entire image is computed once with only 5 x 5 neighborhoods and stored in an array with the size of the image. Computing using a small neighborhood
Figure BDA0001697446790000109
And
Figure BDA00016974467900001010
to prevent excessive smoothing of the image's singular structures. RIBNLM filters having lowerThe complexity is calculated and the blobs, which are structurally similar but have different orientations relative to the reference blob, are given appropriate weights.
Step 5) inverse DNST conversion is carried out on the processed coefficients
The high-frequency sub-band and the low-frequency sub-band after threshold shrinkage in the above steps are subjected to inverse DNST transformation, so that a denoised CT image which is beneficial to analysis of a doctor can be obtained, and the superiority of the algorithm of the invention in denoising the medical CT image is verified by comparing experimental data. The specific algorithm process of the DNST inverse process is described as follows:
s1 inputting the shear coefficient shearlet coeffs (i) epsilon R after DNST processingX*Y*nth；
S2 setting frec∈RX*YRepresenting the reconstructed image sequence;
s3, calculating each index i to be 0, nth]Reconstructed image sequence frequency spectrum f of lower shearlet coeffs (i)recAnd sum frecAccording to the convolution theory and the frame theory
Figure BDA00016974467900001011
S4, inverse transform of DNST is carried out to obtain a reconstructed image sequence frec:＝IFFT(frec)；
The overall steps of the present invention are shown in flow chart form in fig. 5.
Case analysis
The invention researches a discrete inseparable shear wave system on the basis of discrete shear wave transformation by taking a specific medical CT image as an object, adopts a self-adaptive threshold and bivariate shrinkage algorithm in a high-frequency sub-band and adopts rotation invariant bivariate local mean filtering processing in a low-frequency sub-band, and simultaneously shows the superior performance of the invention by comparing with the prior art.
The present invention uses peak signal-to-noise ratio (PSNR) to represent the quality of the reconstructed image, which is defined as follows:
Figure BDA0001697446790000111
where N represents the number of pixels in the image, | | g | | luminanceFRepresenting the frobenius norm, 255 is the maximum value that a pixel can obtain in a grayscale image. The larger the PSNR value is, the better the denoising effect is.
The experiment of the invention adopts medical CT noise image as input data to carry out effective contrast experiment. The decomposition scale is 4, the shear wave level is [1,1,2,2], and the algorithm decomposes the noise map into 48 high frequency subband images and 1 low frequency subband image, which are equal in size to the original image. And (3) carrying out self-adaptive threshold and bivariate shrinkage algorithm on the shear wave coefficient of the high-frequency part, adopting rotation invariant bilateral non-local mean filtering processing on the low-frequency part in the CT according to the characteristics of the low-frequency part, and carrying out shear wave inverse transformation reconstruction on the processed high-frequency and low-frequency data. The experimental results are shown in fig. 6, which compares the experiments with SWT (wavelet transform), FDCT (fast discrete wavelet transform), NSST (non-subsampled shear transform), and FFST (fast finite shear transform). And the denoising effect of the local area in the image is selected for comparison, and the specific details are shown in FIG. 7, so that the superiority of the method can be obviously observed.
As can be seen from table 1, the requirement for the denoising algorithm is higher as the noise variance is larger. On the same noise variance, the algorithm of the invention numerically leads NSST, and has obvious effect more than SWT, FDCT and FFST. And the algorithm of the present invention has more clear detailed features in fig. 6 and 7.
Table 1: PSNR/dB value of different noise reduction algorithms of medical CT image
Algorithm σ＝10 σ＝20 σ＝30 σ＝40 σ＝50
Algorithm of the invention 32.10 28.56 26.91 25.85 24.96
NSST 31.52 28.18 26.40 25.56 24.33
SWT 30.32 26.85 25.24 24.02 23.18
FDCT 29.72 26.84 25.27 24.28 23.58
FFST 29.82 26.62 25.32 24.41 23.76
While the invention has been described in connection with specific embodiments thereof, it will be understood that these should not be construed as limiting the scope of the invention, which is defined in the following claims, and any variations which fall within the scope of the claims are intended to be embraced thereby.

Claims (1)

1. The DNST-based medical CT image denoising method comprises the following steps:
step 1) establishing a medical CT image model; analyzing the noise source in CT, counting abstract noise, and establishing a medical CT model, which specifically comprises the following steps:
statistical noise in CT images, which may also be referred to as quantum noise; this is gaussian additive noise in terms of a mathematical model, which is as follows:
Inoisy(x,y)＝uorignal(x,y)+ηnoisy(x,y) (1)
where (x, y) denotes a pixel point in the image, InosiyRepresenting CT images containing noise, i.e. observed CT images, uorignalRepresenting noiseless CT images, etanoisyRepresenting noise, the distribution of which conforms to a gaussian distribution;
step 2) decomposing the CT image by dispersing the inseparable shear wave; and calculating the shear wave coefficient of the decomposed sub-band, specifically comprising:
s1: inputting a two-dimensional CT image signal f ∈ RX*YA scale parameter J ∈ N, a direction parameter k ∈ NJAnd selecting a directional Filter Direction Filter and a low-pass Filter quad middle Filter;
s2: converting the frequency spectrum f of a CT imagefreq＝FFT(f)；
S3: calculating sub-band i epsilon [0, nth ] under different scales and different directions]Shear wave coefficient under (i) sheearlettCoeffs (i) epsilon RX*Y*nthAccording to the convolution theory of the shear wave domain, the following can be obtained:
Figure FDA0003104454880000011
s4: calculating to obtain discrete unseparatable shear wave coefficients shearlet coeffs (i);
where nth in step S3 represents the redundancy of the entire tightly supported DNST system, which is calculated as follows:
nth＝2*((2*2k[0]+1))+2*((2*2k[1]+1))+...+2*((2*2k[J]+1) (3)
decomposition of medical CT images into f in the frequency domain by DNST decomposition1,f2,...,fnth-1High-frequency CT image and low-frequency CT image f with equal sheet sizenth；
Step 3) on the high frequency shear wave coefficient f1,f2,...,fnth-1Carrying out bivariate shrinkage treatment; through the interrelation of the coefficients among the shear wave coefficient models, the processing is carried out through a bivariate contraction function, and the method specifically comprises the following steps:
utilizing a non-Gaussian bivariate function to count a shear wave coefficient model; in general, for a given original noise-free image, if it is disturbed by additive white gaussian noise, the degraded image is represented in the shear domain as follows:
g＝f+ε (4)
wherein g, f, epsilon respectively represent the observed shear wave coefficient, the original noiseless shear wave coefficient and the noise coefficient; the purpose of denoising is to obtain an estimate of the original coefficients
Figure FDA0003104454880000021
So that the estimated value
Figure FDA0003104454880000022
As close as possible to the originalStarting a noise-free coefficient f; in the proposed method, the maximum a posteriori probability MAP is used to estimate the de-noising coefficients
Figure FDA0003104454880000023
Therefore:
Figure FDA0003104454880000024
bayesian rules can be used to write:
Figure FDA0003104454880000025
the following equation can thus be derived:
Figure FDA0003104454880000026
to describe the dependency between CS and FS in the framework domain, let f1Denotes f2CS, f1Is and f2The frame coefficient at the same spatial position, therefore, f ═ f (f)1,f2)，g＝(g1,g2)，ε＝(ε1,ε2)，f1And f2Are not relevant, but not independent; thus g ═ f + epsilon can be written in the form:
Figure FDA0003104454880000027
assume that the mean and variance of the additive noise are 0 and
Figure FDA0003104454880000028
from pεThe noise probability density function represented by (ε) is as follows:
Figure FDA0003104454880000029
from the distribution of the shear wave coefficients, it is possible to pass through the bivariate probability density function p from a symmetrical circular probability distributionf(f) Fitting the model; thus, it can be expressed as follows:
Figure FDA00031044548800000210
substituting equations (9) and (10) into equation (7) yields:
Figure FDA00031044548800000211
if p isf(f) Being convex, the solution process equation (11) can be converted to the following solution equation:
Figure FDA0003104454880000031
wherein
Figure FDA0003104454880000032
Substituting equation (13) into equation (12), f1The MAP estimator of (a) is a bivariate contraction function:
Figure FDA0003104454880000033
from this equation, it can be seen that for each shear coefficient, there is a corresponding threshold that depends not only on the CS coefficient but also on the FS coefficient;
step 4) on the low frequency shear wave coefficient fnthCarrying out rotation invariant bilateral non-local mean filtering processing; the method specifically comprises the following steps:
using a rotation invariant bilateral non-local mean filter RIBNLM as a post-processing step that optimizes the trade-off between the required smoothness of the individual regions and preserves the details of the image; the proposed filter uses a rotation invariant similarity measure based on intensity values and corresponding neighborhoods to supplement the mean value as
Figure FDA0003104454880000034
And
Figure FDA0003104454880000035
as follows:
Figure FDA0003104454880000036
Figure FDA0003104454880000037
wherein SW represents a search window of size 11 × 11 pixels, wRIBNLM(i, j) is a weight, and the condition 0. ltoreq. w is satisfiedRIBNLM(i, j) is less than or equal to 1 and sigma-wRIBNLM(i,j)＝1，
Figure FDA0003104454880000038
And
Figure FDA0003104454880000039
is the mean of the reference blobs and the blob being processed in the search window SW, and the parameter r is a smoothing parameter; local mean
Figure FDA00031044548800000310
And
Figure FDA00031044548800000311
calculating the whole image once by using 5 multiplied by 5 neighborhoods, and storing the whole image in an array with the image size; computing using a small neighborhood
Figure FDA00031044548800000312
And
Figure FDA00031044548800000313
to prevent excessive smoothing of the singular structure of the image; the RIBNLM filter has a low computational complexity and gives appropriate weights to the blobs, which are structurally similar but have different orientations with respect to the reference blob;
step 5) after filtering the high-frequency shear wave coefficient subjected to bivariate shrinkage and the rotation invariant bivariate non-local mean value subjected to low-frequency processing, carrying out inverse DNST (deep discrete wavelet transform) on the processed coefficient, wherein the specific steps are as follows:
carrying out inverse DNST (digital noise test) transformation on the high-frequency sub-band and the low-frequency sub-band after threshold contraction in the above step to obtain a denoised CT image which is beneficial to analysis of a doctor; the specific algorithm process of the DNST inverse process is described as follows:
t1 input shear coefficient shearlet coeffs (i) epsilon R after DNST processingX*Y*nth；
T2 setting frec∈RX*YRepresenting the reconstructed image sequence;
t3 calculating each index i ∈ [0, nth >]Reconstructed image sequence frequency spectrum f of lower shearlet coeffs (i)recAnd sum frecAccording to the convolution theory and the frame theory
Figure FDA0003104454880000041
T4 obtaining a reconstructed image sequence f by inverse DNST transformationrec:＝IFFT(frec)。
CN201810618315.0A 2018-06-15 2018-06-15 DNST (deep depth transform) -based medical CT (computed tomography) image denoising method Active CN109035156B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810618315.0A CN109035156B (en) 2018-06-15 2018-06-15 DNST (deep depth transform) -based medical CT (computed tomography) image denoising method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810618315.0A CN109035156B (en) 2018-06-15 2018-06-15 DNST (deep depth transform) -based medical CT (computed tomography) image denoising method

Publications (2)

Publication Number Publication Date
CN109035156A CN109035156A (en) 2018-12-18
CN109035156B true CN109035156B (en) 2021-07-30

Family

ID=64609331

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810618315.0A Active CN109035156B (en) 2018-06-15 2018-06-15 DNST (deep depth transform) -based medical CT (computed tomography) image denoising method

Country Status (1)

Country Link
CN (1) CN109035156B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110264429B (en) * 2019-03-25 2023-06-06 山东大学 Image enhancement method based on sparsity and similarity priori
CN110827223B (en) * 2019-11-05 2021-08-24 郑州轻工业学院 CS high-noise astronomical image denoising and reconstructing method combined with fractional order total variation

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103077508A (en) * 2013-01-25 2013-05-01 西安电子科技大学 Transform domain non local and minimum mean square error-based SAR (Synthetic Aperture Radar) image denoising method
CN103093441A (en) * 2013-01-17 2013-05-08 西安电子科技大学 Image denoising method based on non-local means and double variant model of transform domain
US8823374B2 (en) * 2009-05-27 2014-09-02 Siemens Aktiengesellschaft System for accelerated MR image reconstruction
CN106157261A (en) * 2016-06-23 2016-11-23 浙江工业大学之江学院 The shearler of translation invariance converts Medical Image Denoising method

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8823374B2 (en) * 2009-05-27 2014-09-02 Siemens Aktiengesellschaft System for accelerated MR image reconstruction
CN103093441A (en) * 2013-01-17 2013-05-08 西安电子科技大学 Image denoising method based on non-local means and double variant model of transform domain
CN103077508A (en) * 2013-01-25 2013-05-01 西安电子科技大学 Transform domain non local and minimum mean square error-based SAR (Synthetic Aperture Radar) image denoising method
CN106157261A (en) * 2016-06-23 2016-11-23 浙江工业大学之江学院 The shearler of translation invariance converts Medical Image Denoising method

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
《AR Image De-noising Using Local Properties Analysis And Discrete Non-separable Shearlet Transform》;Ye Yuan, Liangzhuo Xie, Yewen Zhu, Sheng Wang, Zhemin Zhuang;《2017 10th International Congress on Image and Signal Processing, BioMedical Engineering and Informatics (CISP-BMEI 2017)》;20171231;全文 *
《Nonseparable Shearlet Transform》;Wang-Q Lim;《IEEE Trans Image Process，2013，22( 5) : 2056-2065》;20130531;全文 *
《基于Shearlet变换的图像融合与去噪方法研究》;高国荣;《中国博士学位论文全文数据库-信息科技辑》;20160331;全文 *
《基于剪切波变换的图像处理技术的研究》;徐荣;《中国优秀硕士学位论文全文数据库-信息科技辑》;20161130;全文 *

Also Published As

Publication number Publication date
CN109035156A (en) 2018-12-18

Similar Documents

Publication Publication Date Title
Sagheer et al. A review on medical image denoising algorithms
Diwakar et al. A review on CT image noise and its denoising
Andria et al. Linear filtering of 2-D wavelet coefficients for denoising ultrasound medical images
Wu et al. Evaluation of various speckle reduction filters on medical ultrasound images
CN109598680B (en) Shear wave transformation medical CT image denoising method based on rapid non-local mean value and TV-L1 model
CN109961411B (en) Non-downsampling shear wave transformation medical CT image denoising method
Kaur et al. A comprehensive review of denoising techniques for abdominal CT images
Sudha et al. Speckle noise reduction in ultrasound images using context-based adaptive wavelet thresholding
Mishro et al. A survey on state-of-the-art denoising techniques for brain magnetic resonance images
EP2567357A1 (en) A method and device for estimating noise in a reconstructed image
CN109003232B (en) Medical MRI image denoising method based on frequency domain scale smoothing Shearlet
Farouj et al. Hyperbolic Wavelet-Fisz denoising for a model arising in Ultrasound Imaging
CN109035156B (en) DNST (deep depth transform) -based medical CT (computed tomography) image denoising method
Diwakar et al. CT Image noise reduction based on adaptive wiener filtering with wavelet packet thresholding
CN109559283B (en) DNST domain bivariate shrinkage and bilateral non-local mean filtering-based medical PET image denoising method
Dawood et al. The importance of contrast enhancement in medical images analysis and diagnosis
Singh et al. A new local structural similarity fusion-based thresholding method for homomorphic ultrasound image despeckling in NSCT domain
CN108846813A (en) The medicine CT image denoising method of frame and NSST is decomposed based on MFDF
CN109584322B (en) Shearlet medical PET image denoising method based on frequency domain direction smoothing
CN112734663A (en) Profile wave transformation medical CT image noise reduction method based on self-adaptive threshold
Devi et al. An improved adaptive wavelet shrinkage for ultrasound despeckling
CN109767410B (en) Lung CT and MRI image fusion algorithm
Biradar et al. Echocardiographic image denoising using extreme total variation bilateral filter
Sneha et al. Performance analysis of wiener filter in restoration of COVID-19 chest X-ray images, ultrasound images and mammograms
CN102402783A (en) Method for processing spots of three-dimensional ultrasonic image

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant