WO2013014110A1 - Procédé de réduction du bruit dans une séquence d'images fluoroscopiques par filtrage temporel et spatial - Google Patents

Procédé de réduction du bruit dans une séquence d'images fluoroscopiques par filtrage temporel et spatial Download PDF

Info

Publication number
WO2013014110A1
WO2013014110A1 PCT/EP2012/064347 EP2012064347W WO2013014110A1 WO 2013014110 A1 WO2013014110 A1 WO 2013014110A1 EP 2012064347 W EP2012064347 W EP 2012064347W WO 2013014110 A1 WO2013014110 A1 WO 2013014110A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
coefficients
coefficient
value
acquired
Prior art date
Application number
PCT/EP2012/064347
Other languages
English (en)
Inventor
Jérémie Pescatore
Carole AMIOT
Original Assignee
Thales
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 Thales filed Critical Thales
Priority to US14/234,292 priority Critical patent/US9058654B2/en
Priority to EP12737811.5A priority patent/EP2734975B1/fr
Priority to JP2014520683A priority patent/JP6120843B2/ja
Publication of WO2013014110A1 publication Critical patent/WO2013014110A1/fr

Links

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/10116—X-ray image
    • G06T2207/10121—Fluoroscopy
    • 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/20004—Adaptive image processing
    • G06T2207/20012—Locally adaptive
    • 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/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
    • 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/20172—Image enhancement details
    • G06T2207/20182—Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing

Definitions

  • the invention is in the field of X-ray imaging and, more specifically, that of medical fluoroscopy imaging. It relates to a method of reducing noise in a sequence of fluoroscopic images acquired by an X-ray detector.
  • X-ray fluoroscopy imaging makes it possible to provide a flow of images from a patient to a physician during so-called minimally invasive surgical procedures, for example the chemoembolization of the liver. vertebroplasty, catheterization of aneurysms or treatment of vascular stenosis.
  • the images can notably help guide surgical instruments.
  • a fluoroscopically guided intervention is typically performed by passing a catheter within the vascular network of the patient.
  • a contrast agent may or may not have been previously injected into the vascular network in order to opacify it and improve the visibility of the vessels.
  • Such an intervention is generally relatively long and the dose of X-rays to which the patient is subjected must be limited in order to avoid causing tissue lesions or burns. Because of this limitation of the X-ray dose used, the fluoroscopic images obtained have a relatively high noise level, and therefore a relatively low contrast-to-noise ratio, making the images difficult to read.
  • Filtering treatments are implemented in order to reduce the quantum noise present in these images and increase their contrast-to-noise ratio.
  • FNR treatment A treatment known under the name of fluoroscopic noise reduction treatment, or FNR treatment according to the abbreviation of the English expression "Fluoroscopy Noise Reduction”, is thus generally performed on the images acquired by the X-ray detector.
  • the objective of this FNR treatment is to filter the noise present in the image while preserving the contrast of the information present in this image.
  • the FNR is performed by applying a temporal filter to the regions of the images where there is no displacement. The existence or absence of displacements in regions of the image is detected from the individual variation of the intensity of each pixel considered separately.
  • a pixel is considered to be in motion when its intensity variation between two images exceeds a threshold linked to the standard deviation of the noise.
  • a recursive temporal filtering treatment is applied to the pixels determined to be fixed.
  • the FNR treatment has the advantage of not coloring the noise spatially. In other words, it does not reveal fictitious objects resulting from groupings of pixels of intensities close to each other.
  • the FNR treatment has a relatively limited denoising capacity and tends to remove objects of interest in the image.
  • spatial filtering has interesting properties.
  • filtering by image transformation in the wavelet domain is one of the most used techniques.
  • the wavelet theory introduced in the 1980s, was developed from the Fourier transform.
  • the Fourier transform is used to describe many physical phenomena. It is indeed a powerful mathematical tool to describe many physical phenomena, including physical phenomena modeled by stationary signals.
  • the representations of physical phenomena modeled by non-stationary signals are not generally satisfactory.
  • the Fourier transform is unable to locate the portions of the signal where the frequency suddenly changes.
  • a Fourier transform with windows has been developed. It consists of multiplying the signal by a window of given dimension.
  • the analysis thus becomes local, the size of the window determining the temporal resolution obtained, or the spatial resolution in the case of two-dimensional signals.
  • the wavelet theory is based on the Fourier transform with windows. It differs in that it allows, during the transformation, to change the size of the window. For this reason, we speak of multiresolution analysis or multiscale analysis. Another difference is that the signal is no longer broken down into a sum sine and cosine, but according to functions called "wavelets".
  • the wavelets have the advantage of transforming the image into a so-called hollow representation.
  • hollow representation is meant a representation requiring a small number of parameters to faithfully represent the image.
  • An object of the invention is notably to provide a denoising treatment adapted to fluoroscopic images whose contrast to noise ratio is relatively low.
  • the denoising treatment must have a high denoising capacity, generate a low calculation cost, limit the loss of objects of interest in the image, and not introduce artifacts.
  • the subject of the invention is a method for reducing noise in a sequence of fluoroscopic images acquired by an X-ray detector, each image being formed in the spatial domain of a matrix of pixels, each pixel having a value representative of a signal level, the method comprising the following successive steps for each image:
  • ⁇ apply temporal filtering on the image acquired at an instant n, the application of said temporal filtering comprising the following substeps for each pixel of the image:
  • ⁇ apply spatial filtering on the image acquired at time n, the application of said spatial filtering comprising the following substeps:
  • thresholding the coefficients of the image by a thresholding function, the thresholding function canceling the coefficients lower than a third predetermined threshold, and keeping or adjusting the coefficients greater than the third predetermined threshold,
  • the stage of transformation of the image acquired at time n in the field of curvelets advantageously uses a transform into curvelets with 6 or 9 scales, and / or 16 orientations.
  • the step of transforming the image acquired at time n in the curvelet domain uses a discrete curvature transform by USFFT.
  • the step of transforming the image acquired at time n in the curvelet domain uses a discrete curvelet transformation by wrapping.
  • the thresholding function of the spatial filtering may be a hard thresholding function, that is to say a thresholding function according to which the coefficients are canceled if they are lower than the third predetermined threshold, and otherwise preserved.
  • the third predetermined threshold T s of the spatial filtering can be determined by the following successive steps:
  • ⁇ determine one of the coefficients c i Y k by the relation: where M and P designate the dimensions of the image Y and / denotes the orientations of the curvelets,
  • ⁇ is the standard deviation of the noise in the image acquired at time n and on which the temporal filtering has been applied
  • the factor r takes the value 0 for the coarsest scale of the curvelets, the value 4 for the finest scale, and 1 for the other scales.
  • the step of applying the temporal filtering on the image acquired at the instant n comprises, following the substep of determining the difference between the value of the pixel considered in the image acquired at time n and the value of the corresponding pixel in the image acquired at time n-1 on which temporal filtering has been applied, the following additional sub-step:
  • the noise reduction method according to the invention successively comprises temporal filtering and spatial filtering by curvelet transformation.
  • the order in which the filtering is carried out is of importance with regard to the coloration of the noise.
  • spatial filtering by curvelet transformation introduces artifacts if it is not preceded by temporal filtering.
  • the invention has the particular advantage of not requiring a priori forms for spatial filtering.
  • FIG. 1 represents an example of points taken into account for determining the contrast of an object in a fluoroscopic image
  • FIG. 2 represents possible steps of the noise reduction method according to the invention
  • FIGS. 3A and 3B show, by curves, a hard thresholding function and a soft thresholding function, respectively;
  • FIG. 4 illustrates, by a graph, the performance of this noise reduction method for a sequence of fluoroscopic images representing an object having a translational and rotational movement
  • FIG. 5 represents an example of sub-steps of the noise reduction method making it possible to perform spatial filtering
  • FIG. 6 represents an example of a thresholding function used for the spatial filtering of FIG. 5.
  • CNR contrast-to-noise ratio
  • FIG. 1 represents the points considered (in white) for determining the contrast of a stent, better known by the English name of "stent".
  • Figure 2 shows possible steps for the noise reduction method in the fluoroscopic image sequence. The method is described for the image x n acquired at time n. However, it can be applied successively to each of the images in the sequence.
  • the image x n is processed by a preprocessing chain.
  • Pretreatment may include groupings of adjacent pixels in the image.
  • the pixels can be grouped for example according to the desired resolution or the desired noise level.
  • the pretreatment may also comprise a phase of calibration of the X-ray detector. In particular, the gain and the offset to be applied to each pixel of the detector may be determined.
  • the calibration may also include an identification of the pixels taking a random value or a constant value regardless of the received radiation. These defective pixels can be corrected by interpolating their value by the values of the neighboring pixels.
  • the pretreatment may also include rescaling the images, i.e., an overall adjustment of the value of all pixels of each image, so that the average pixel value is constant for all footage of the sequence.
  • the preprocessing may comprise a substep of adaptation of the images.
  • the noise present in each image x n is a fish noise. Such noise is unsuitable for noise reduction treatments.
  • each image x n is a representative image of X-ray attenuations, while an image of radiological thicknesses is more natural to analyze. Ascombe's formula makes it possible to obtain an image of the radiological thicknesses with a pseudo-Gaussian noise:
  • x 0 n (m, p) is the value of the pixel (m, p) before the substep of adaptation of the image x n .
  • the image x n is filtered by a temporal filter.
  • the temporal filtering is for example a filtering called "Fluoroscopic Noise Reduction" (FNR). This filtering is described in particular in US Pat. No. 6,314,160 B1.
  • FNR Fluoroscopic Noise Reduction
  • a (m, p) the difference between the value x n (m, p) in the image x n acquired at time n and the value y n _i (m, p) in filtered image n. ! corresponding to the acquisition at time n-1:
  • FNR filtering is based on two modes, called adaptive mode and non-adaptive mode.
  • the adaptive mode is called because of its consideration of the movement of objects.
  • FNR filtering is to compare pixel by pixel, the difference ⁇ at predetermined thresholds.
  • the first threshold Si corresponds to a relatively small difference ⁇ .
  • the second threshold S 2 corresponds to a relatively high difference ⁇ .
  • FNR filtering is then parameterized by coefficients a and ⁇ as a function of the result of this difference ⁇ and the mode considered.
  • the strength of the filter is:
  • ⁇ DD ⁇ (1 - a)
  • the thresholds Si and S 2 are defined as follows:
  • the thresholds Si and S 2 are defined as follows:
  • FNR filtering itself consists in determining a new value y n for each pixel of the image from its value x n and the coefficients a and ⁇ :
  • y n a. x n + (1 - a) n _ ! + ⁇
  • the filtering step FNR is shown for the adaptive mode.
  • the absolute value of the difference ⁇ is compared for each pixel x n at the threshold S lt, namely the standard deviation o. If the absolute value of the difference ⁇ is less than the S lt threshold the pixel is probably a high noise level and said filter is applied in a sub-step 222 with the corresponding parameters. Otherwise, in a sub-step 223, the absolute value of the difference ⁇ is compared for each pixel x n at the threshold S 2 , namely 2o.
  • the pixel considered corresponding probably to a point of the moving object and so-called weak filtering is applied in a substep 224 with the corresponding parameters.
  • the absolute value of the difference ⁇ is between the thresholds 5 X and S 2 , an uncertainty exists.
  • an average filtering is applied in a sub-step 225 with the corresponding parameters.
  • a third step 23 the image y n obtained after time filtering is filtered by a spatial filter based on the curvelet transform.
  • the theory of curvelets uses the principle of multiresolution or multiscale spaces, a scale corresponding to a partition level of space.
  • the space is partitioned into curvelets whose envelopes are of dimensions 2 j x 2 d / 2 , where j is a positive integer denoting the scale.
  • Each curvelet depends on the number of scales considered, its position in the image and its orientation. For each scale, the best possible approximation of the image is calculated. From one scale to another, only the improvements obtained are represented. The number of scales determines the quality of the reconstruction of the image.
  • the image y n is transformed from the spatial domain to the curvelet domain by a transformation into curvelets.
  • a discrete curvelet transformation is used.
  • Such a transformation is notably described in E. Candès, Demanet L, Donoho D. and Ying L, "Fast Discrete Curvelet Transform", Multiscale Model. Simul., Vol. 5, no. 3, 2006, pp. 861 - 899.
  • the curvelets can not be sampled according to the Fourier space grids.
  • the curvelet transform of the image y n can be represented by a set of coefficients c (j, l, k), where the exponent D refers to the discrete transform, and where j, l and k denote respectively the scale, the orientation and the position of the curvelet considered.
  • the number of scales used is equal to 9.
  • the number of scales used is equal to 6. This last embodiment implies a smaller number of coefficients to calculate and therefore a lower processing time.
  • the number of orientations used is advantageously equal to 16. This number presents a good compromise between the resolution in orientation and the computational cost of the curvelets, in particular for real images of which curves require a large number of orientations to approach them. Moreover, the number 16 being a multiple of 4, it facilitates the transformation into curvelets by wrapping.
  • the coefficients c (j, l, k) of the image y n are thresholded by a thresholding function.
  • the thresholding step 232 is based on the assumption that some coefficients c (j, l, k) represent the noise while others contain the useful information.
  • the thresholding step 232 therefore consists of identifying and deleting the coefficients that do not contain useful information.
  • the thresholding function is characterized on the one hand by its shape and on the other hand by its threshold T. There are mainly two thresholding functions, the so-called soft thresholding function and the so-called hard thresholding function. By simplifying the notation of the coefficients by c n and calling c n the thresholded coefficients, the hard thresholding function can be modeled as follows:
  • the soft thresholding function can be modeled as follows: if not
  • FIG. 3A shows a curve 31 the hard thresholding function and FIG. 3B shows a curve 32 the soft thresholding function.
  • the soft thresholding function has a greater denoising capacity than the hard thresholding function.
  • it modifies the value of the pixels according to the threshold T chosen. This results in a loss of precision on certain details of the image.
  • the determination of the threshold T is a crucial step because it determines the performance of spatial filtering, both in terms of denoising capacity and the introduction of artifacts.
  • the threshold T is determined as being the universal threshold ⁇ ⁇ described in particular in DL Donoho and Johnstone, IM, "Ideal spatial adaptation via wavelet shrinkage", Biometrica, vol. 81, 1994, pp. 425-455.
  • This threshold ⁇ ⁇ is defined by the following relation:
  • N is the number of samples, i.e. here the number of pixels in the image y n .
  • the universal threshold ⁇ ⁇ is not optimal and can introduce artifacts. However, thresholds proportional to r universal threshold there may be good results if the proportionality factor is chosen.
  • the threshold T is determined by the method described in JL Starck, Candès IJ and Donoho L, "The Curve Transform for Image Denoising", IEEE Transactions on Image Processing, vol. 1 1, no. 6, 2002, pp. 670-684.
  • This method consists in creating an image, denoted by Y, of the same dimensions as the image y n to be denoised and of uniform value 1 on all the pixels.
  • the Y image is transformed by the Fourier transform.
  • the coordinate frequency (0,0) in the image Y is centered and its amplitude is normalized by the dimensions of the image.
  • F be the amplitude of the central frequency. Its normalized amplitude, noted
  • M and P designate the dimensions of the image Y, that is to say also of the image x n and of the image y n .
  • T s The threshold of this method, denoted T s, depends on the scale j and on the orientation l. It is determined according to the standard of the coefficients c i Y k of the curvelets of the image Y for a scale j and a direction l data. This standard c i Y is determined by the following relation:
  • the thresholding function is a hard thresholding function, defined as follows:
  • the image y n is transformed from the domain of curvelets to the spatial domain (matrix) by a transformation into inverse curvelets.
  • the threshold coefficients are used.
  • the image obtained after temporal filtering and spatial filtering is denoted z n .
  • FIG. 4 graphically illustrates the performance of the noise reduction method according to the invention for a sequence of images representing a stent having translational and rotational movement.
  • the x-axis represents the average CNR of the x n sequence of images and the y-axis represents the average CNR of the image sequence obtained after noise reduction processing.
  • a first curve 41 represents the average CNR of the image sequence z n obtained by the noise reduction method according to the invention, that is to say by a temporal filtering followed by a spatial filtering by transformation into curvelets.
  • a second curve 42 represents the average CNR of the same sequence of images after processing by the same temporal filtering followed by spatial filtering by wavelet transformation.
  • the graph shows that for a sequence of images to be processed with a relatively low CNR (0.5-0.7-1), the noise reduction method according to the invention shows a significant improvement in terms of CNR.
  • the noise reduction is even better for a sequence of images with a higher CNR (1, 5 - 2 - 2.5).
  • the noise reduction method according to the invention makes it possible to significantly reduce the introduction of artifacts.
  • the spatial filtering step 23 described above makes it possible to effectively reduce the noise level in the images.
  • this step has the disadvantage of introducing artifacts appearing in the form of curvelets, that is to say objects with elongated and undulating shapes. Indeed, these artifacts are due to coefficients having a large amplitude, but representing only noise. They are present only when the noise level is high, which is the case of fluoroscopic images.
  • One solution is to set relatively low thresholds, that reduces artifacts but limits performance in terms of noise reduction.
  • the following description describes another method for performing the spatial filtering step to significantly reduce the presence of artifacts without reducing the performance of the filtering.
  • This method is based on the correlation that exists between a given coefficient and the coefficients that are related to it, namely its parent coefficient, its neighboring coefficients and its cousin coefficients.
  • a coefficient represents the useful signal
  • the coefficients that are bound to it have a high probability of also representing the useful signal.
  • the modulus of a coefficient is large because of the presence of an artifact, it is unlikely that the coefficients related to it are also affected by this artifact.
  • FIG. 5 represents another example of sub-steps for performing the spatial filtering step 23.
  • the image y n obtained after temporal filtering is transformed from the spatial domain to the curvelet domain.
  • Transform curvelets image y n is represented by rated ° coefficients c (j, l, k), or c n in simplified form.
  • the transformation into curvelets is done with 6 scales and 12 orientations.
  • T x , T 2 and T 3 are determined such that T l ⁇ T 2 ⁇ T 3 .
  • These thresholds are determined according to the standard deviation a sfl of the subband of the coefficients.
  • the subband of a given coefficient c ° (j, l, k) groups all the coefficients having the same scale j and the same orientation /. It is therefore the set of positions for the scale and orientation considered.
  • the thresholds ⁇ ⁇ , T 2 and T 3 define a first zone for the values below the threshold ⁇ ⁇ , a second zone for the values between the thresholds T x and T 2 , a third zone for the values between the thresholds T 2 and T 3 , and a fourth zone for values above the threshold T 3 .
  • the coefficients are corrected according to their parent coefficient and their neighboring coefficients. Of preferably, only the uncorrelated neighbor coefficients are considered.
  • c ° (j, l, k), or c n in simplified form, the corrected coefficients.
  • the second case is where the modulus of the parent coefficient c ° (jl, l, k) and the mean of the modules of the neighboring coefficients c ° (j, l, k ') lie in the same area, different from that of the module the coefficient considered c (/, /, &).
  • the modulus of the corrected coefficient c ° (j, l, k) then takes the minimum value of the area in which the modulus of the parent coefficient is situated and the average of the modules of the neighboring coefficients if the modulus of the coefficient considered c ° (j, l, k) is in a zone with lower values, and the maximum value of this zone otherwise.
  • the third case is the one where the modulus of the coefficient considered c ° (j, l, k), the modulus of its parent coefficient c ° (jl, l, k), and the mean of the modules of its neighboring coefficients c ° (j , l, k ') are in three distinct zones.
  • the modulus of the corrected coefficient c (j, l, k) then takes the minimum value of the area in which the average of the modules of the neighboring coefficients is situated if the modulus of the considered coefficient c ° (j, l, k) lies in a zone with lower values, and the maximum value of this zone otherwise.
  • the third case favors the zone of the uncorrelated neighboring coefficients to the detriment of the zone of the parent coefficient because at the coarsest scale, the coefficients do not have a parent. It would nevertheless be possible to correct the coefficients of other scales according to the area in which the module of their parent coefficient is located.
  • the corrected coefficients c n are thresholded by a thresholding function.
  • a thresholding function We write c ° (j, l, k), or c n in simplified form, the thresholded corrected coefficients.
  • FIG. 6 represents the thresholding function corresponding to the threshold values T x , T 2 and T 3 considered. This function is modeled as follows:
  • f ⁇ c n is a sigmoid defined as follows:
  • the parameters m and ⁇ are determined so that the thresholding function is continuous. They are thus defined as follows:
  • the thresholding function thus cancels the coefficients whose modulus is included in the first zone 61, it applies a soft thresholding to the coefficients whose module is included in the second zone 62, it applies a hard thresholding to the coefficients whose module is included in the fourth zone 64, and it applies a transition threshold to the coefficients whose module is included in the third zone 63.
  • the soft thresholding makes it possible to reduce the noise level, but it spreads the outlines. Hard thresholding keeps the contours well but is less efficient in terms of noise reduction. The joint use of these two types of thresholding makes it possible to combine their respective advantages according to the zones.
  • threshold T 4 is defined from the standard deviation a sfl of the subband of the coefficient c "considered:
  • a fifth substep 55 the image y n is transformed from the curvelet domain to the spatial domain by a transformation into inverse curvelets.
  • the corrected threshold coefficients are used.
  • the image obtained after temporal filtering and spatial filtering is denoted z ' n .
  • the spatial filtering described with reference to Figure 5 uses both the parent coefficients and the uncorrelated neighbor coefficients. Spatial filtering could also use correlated neighboring coefficients, or the set of neighboring coefficients. It could also use only the parent coefficients or the neighboring coefficients, or combine these coefficients with the cousin coefficients. These coefficients are those having the same scale and the same position as the reference coefficient.
  • the spatial filtering could also not include a sub-step of correction of the coefficients, the sub-step 54 of thresholding being in this case carried out directly on the coefficients c n of the transforms in curvelets.

Landscapes

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

Abstract

L'invention se situe dans le domaine de l'imagerie médicale par fluoroscopie. Elle concerne un procédé de réduction du bruit dans une séquence d'images fluoroscopiques acquises par un détecteur à rayons X. Le procédé selon l'invention comprend les étapes successives suivantes pour chaque image (xn) de la séquence : ■ appliquer un filtrage temporel (22) sur l'image (xn) acquise à un instant n, ledit filtrage temporel comportant un traitement de réduction du bruit fluoroscopique (FNR), ■ appliquer un filtrage spatial (23) sur l'image (yn) acquise à l'instant n, ledit filtrage spatial comportant les sous-étapes suivantes : - transformer (231) l'image (yn) acquise à l'instant n du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image transformée étant représentée par un ensemble de coefficients, - seuiller (232) les coefficients de l'image (yn) par une fonction de seuillage, la fonction de seuillage annulant les coefficients inférieurs à un troisième seuil prédéterminé, et conservant ou ajustant les coefficients supérieurs au troisième seuil prédéterminé, - transformer (233) l'image (zn) dont les coefficients ont été seuillés du domaine des curvelets au domaine spatial par une transformation en curvelets inverse.

Description

PROCÉDÉ DE RÉDUCTION DU BRUIT DANS UNE SÉQUENCE D'IMAGES FLUOROSCOPIQUES PAR FILTRAGE TEMPOREL ET
SPATIAL
L'invention se situe dans le domaine de l'imagerie par rayons X et, plus spécifiquement, celui de l'imagerie médicale par fluoroscopie. Elle concerne un procédé de réduction du bruit dans une séquence d'images fluoroscopiques acquises par un détecteur à rayons X.
Dans le domaine de l'imagerie médicale, l'imagerie par fluoroscopie à rayons X permet de fournir un flux d'images d'un patient à un médecin lors d'interventions chirurgicales dites minimalement invasives, par exemple la chimio-embolisation du foie, la vertébroplastie, la cathérisation des anévrismes ou le traitement des rétrécissements vasculaires. Les images peuvent notamment aider le guidage d'instruments chirurgicaux. Une intervention guidée par fluoroscopie est typiquement réalisée en faisant passer un cathéter à l'intérieur du réseau vasculaire du patient. Un agent de contraste peut ou non avoir été injecté préalablement dans le réseau vasculaire afin de l'opacifier et d'améliorer la visibilité des vaisseaux. Une telle intervention est généralement relativement longue et la dose de rayons X auxquels le patient est soumis doit être limitée afin d'éviter d'occasionner des lésions ou des brûlures des tissus. Du fait de cette limitation de la dose de rayons X utilisée, les images fluoroscopiques obtenues comportent un niveau de bruit relativement élevé, et donc un rapport contraste sur bruit relativement faible, rendant les images difficilement lisibles.
Des traitements de filtrage sont mis en œuvre afin de diminuer le bruit quantique présent dans ces images et augmenter leur rapport contraste sur bruit.
Un traitement connu sous le nom de traitement de réduction du bruit fluoroscopique, ou traitement FNR selon le sigle de l'expression anglo- saxonne "Fluoroscopie Noise Réduction", est ainsi généralement réalisé sur les images acquises par le détecteur à rayons X. L'objectif de ce traitement FNR est de filtrer le bruit présent dans l'image tout en préservant le contraste de l'information présente dans cette image. Le traitement FNR est réalisé en appliquant un filtre temporel aux régions des images où il n'y a pas de déplacement. L'existence ou l'absence de déplacements dans des régions de l'image est détectée à partir de la variation individuelle de l'intensité de chacun des pixels considérés séparément. Un pixel est considéré comme étant en déplacement lorsque sa variation d'intensité entre deux images dépasse un seuil lié à l'écart type du bruit. Ces pixels déterminés comme étant en déplacement ne sont pas, ou très peu, filtrés. En revanche, un traitement par filtrage temporel récursif est appliqué aux pixels déterminés comme étant fixes. Le traitement FNR présente l'avantage de ne pas colorer le bruit spatialement. Autrement dit, il ne fait pas apparaître d'objets fictifs résultant de groupements de pixels d'intensités voisines entre elles. Cependant, le traitement FNR présente une capacité de débruitage relativement limitée et tend à faire disparaître les objets d'intérêts dans l'image.
A côté du traitement par filtrage temporel FNR, les traitements par filtrage spatial présente des propriétés intéressantes. Parmi les différentes variantes de traitement par filtrage spatial, le filtrage par transformation de l'image dans le domaine des ondelettes est l'une des techniques les plus utilisées. La théorie des ondelettes, introduite dans les années 1980, a été développée à partir de la transformée de Fourier. La transformée de Fourier est utilisée pour décrire de nombreux phénomènes physiques. Elle constitue en effet un outil mathématique performant pour décrire de nombreux phénomènes physiques, notamment les phénomènes physiques modélisés par des signaux stationnaires. En revanche, les représentations des phénomènes physiques modélisés par des signaux non stationnaires ne sont généralement pas satisfaisantes. En particulier, la transformée de Fourier est incapable de localiser les portions du signal où la fréquence change subitement. Pour y remédier, une transformée de Fourier à fenêtres a été développée. Elle consiste à multiplier le signal par une fenêtre de dimension donnée. L'analyse devient donc locale, la dimension de la fenêtre déterminant la résolution temporelle obtenue, ou la résolution spatiale dans le cas de signaux à deux dimensions. La théorie des ondelettes s'appuie sur la transformée de Fourier à fenêtres. Elle s'en différencie de par le fait qu'elle permet, au cours de la transformation, de changer la dimension de la fenêtre. Pour cette raison, on parle d'analyse multirésolution ou analyse multi-échelle. Une autre différence est que le signal n'est plus décomposé en une somme de sinus et de cosinus, mais selon des fonctions appelées "ondelettes". Les ondelettes présentent l'avantage de transformer l'image en une représentation dite creuse. Par représentation creuse, on entend une représentation nécessitant un faible nombre de paramètres pour représenter fidèlement l'image. Cependant, les ondelettes sont mal adaptées à la représentation de discontinuités de type ligne, telles que les contours d'objets, car leur forme ne leur permet pas d'exploiter la forme des objets. En conséquence, l'information sur ces discontinuités est contenue dans les coefficients de nombreuses échelles. Pour obtenir une faible erreur de reconstruction, il faut donc calculer un grand nombre de coefficients, entraînant un coût de calcul important. Pour une application de fluoroscopie où les images doivent être traitées "en temps réel", c'est-à-dire dans des durées inférieures à la période de rafraîchissement des images, ce coût de calcul implique des ressources matérielles importantes.
D'autres théories ont été développées pour réduire la redondance des informations. Il s'agit notamment des ondelettes isotropes et anisotropes, ou encore des arbres duaux orthogonaux ou non.
Au début des années 2000, la théorie des curvelets est apparue. Cette théorie est également fondée sur le principe de l'analyse multi- échelles. La géométrie des curvelets leur permet de tirer profit de la régularité des contours des objets. Une image peut ainsi être représentée par un nombre relativement faible de coefficients. Cependant, la transformation en curvelets présente l'inconvénient d'introduire des artéfacts dans l'image, notamment une coloration spatiale du bruit. L'introduction d'artéfacts est d'autant plus forte que le rapport contraste sur bruit est faible. Or, dans le domaine de l'imagerie par fluoroscopie, le rapport contraste sur bruit est souvent inférieur à 1 . La présence d'artéfacts apparaît donc comme étant inévitable et a limité jusque là son utilisation à des images dont le rapport contraste sur bruit est bien plus élevé.
Un but de l'invention est notamment de proposer un traitement de débruitage adapté aux images fluoroscopiques dont le rapport contraste sur bruit est relativement faible. En particulier, le traitement de débruitage doit présenter une forte capacité de débruitage, engendrer un faible coût de calcul, limiter la perte des objets d'intérêt dans l'image, et ne pas introduire d'artefacts. A cet effet, l'invention a pour objet un procédé de réduction du bruit dans une séquence d'images fluoroscopiques acquises par un détecteur à rayons X, chaque image étant formée dans le domaine spatial d'une matrice de pixels, chaque pixel ayant une valeur représentative d'un niveau de signal, le procédé comprenant les étapes successives suivantes pour chaque image :
■ appliquer un filtrage temporel sur l'image acquise à un instant n, l'application dudit filtrage temporel comportant les sous-étapes suivantes pour chaque pixel de l'image :
- déterminer une différence entre la valeur du pixel considéré dans l'image acquise à l'instant n et la valeur du pixel correspondant dans l'image acquise à un instant n-1 sur laquelle le filtrage temporel a été appliqué,
- si la différence est inférieure à un premier seuil prédéterminé, corriger la valeur du pixel correspondant avec un premier ensemble de paramètres,
- si la différence est supérieure à un deuxième seuil prédéterminé, corriger la valeur du pixel correspondant avec un deuxième ensemble de paramètres,
■ appliquer un filtrage spatial sur l'image acquise à l'instant n, l'application dudit filtrage spatial comportant les sous-étapes suivantes :
- transformer l'image acquise à l'instant n du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image dans le domaine des curvelets étant représentée par un ensemble de coefficients,
- seuiller les coefficients de l'image par une fonction de seuillage, la fonction de seuillage annulant les coefficients inférieurs à un troisième seuil prédéterminé, et conservant ou ajustant les coefficients supérieurs au troisième seuil prédéterminé,
- transformer l'image dont les coefficients ont été seuillés du domaine des curvelets au domaine spatial par une transformation en curvelets inverse. L'étape de transformation de l'image acquise à l'instant n dans le domaine des curvelets utilise avantageusement une transformée en curvelets à 6 ou 9 échelles, et/ou 16 orientations. Selon une première forme particulière de réalisation, l'étape de transformation de l'image acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets discrète par USFFT.
Selon une deuxième forme particulière de réalisation, l'étape de transformation de l'image acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets discrète par wrapping.
La fonction de seuillage du filtrage spatial peut être une fonction de seuillage dur, c'est-à-dire une fonction de seuillage selon laquelle les coefficients sont annulés s'ils sont inférieurs au troisième seuil prédéterminé, et conservés sinon.
Le troisième seuil prédéterminé Ts du filtrage spatial peut être déterminé par les étapes successives suivantes :
■ créer une image Y de mêmes dimensions que l'image acquise à l'instant n et de valeur uniforme 1 sur tous les pixels,
■ transformer l'image Y par une transformée de Fourier,
■ normaliser l'amplitude F du pic dans la transformée de Fourier de l'image Y par la relation :
6 F
\F\ =— ,
5 x P
où M et P désignent les dimensions de l'image Y,
■ transformer la transformée de Fourier de l'image Y du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image dans le domaine des curvelets étant représentée par un ensemble de coefficients cjl k,
■ déterminer un des coefficients ci Y k par la relation :
Figure imgf000007_0001
où M et P désignent les dimensions de l'image Y et / désigne les orientations des curvelets,
■ déterminer le seuil prédéterminé Ts par la relation :
Ts = r. c] . o
où σ est l'écart type du bruit dans l'image acquise à l'instant n et sur laquelle le filtrage temporel a été appliqué, et où le facteur r prend la valeur 0 pour l'échelle la plus grossière des curvelets, la valeur 4 pour l'échelle la plus fine, et la valeur 1 pour les autres échelles.
Selon une forme particulière de réalisation, l'étape d'application du filtrage temporel sur l'image acquise à l'instant n comprend, à la suite de la sous-étape de détermination de la différence entre la valeur du pixel considéré dans l'image acquise à l'instant n et la valeur du pixel correspondant dans l'image acquise à l'instant n-1 sur laquelle le filtrage temporel a été appliqué, la sous-étape supplémentaire suivante :
■ si la différence est comprise entre le premier seuil prédéterminé et le deuxième seuil prédéterminé, corriger la valeur du pixel correspondant avec un troisième ensemble de paramètres.
En résumé, le procédé de réduction de bruit selon l'invention comporte successivement un filtrage temporel et un filtrage spatial par transformation en curvelets. L'ordre dans lequel sont effectués les filtrages présente une importance en ce qui concerne la coloration du bruit. En effet, le filtrage spatial par transformation en curvelets introduit des artéfacts s'il n'est pas précédé du filtrage temporel. Par ailleurs, l'invention a notamment pour avantage de ne pas nécessiter d'à priori de formes pour le filtrage spatial.
L'invention sera mieux comprise et d'autres avantages apparaîtront à la lecture de la description qui va suivre, faite en regard de dessins annexés sur lesquels :
- la figure 1 représente un exemple de points pris en compte pour déterminer le contraste d'un objet dans une image fluoroscopique ;
- la figure 2 représente des étapes possibles du procédé de réduction de bruit selon l'invention ; - les figures 3A et 3B représentent par des courbes une fonction de seuillage dur et une fonction de seuillage doux, respectivement ;
- la figure 4 illustre, par un graphique, la performance de ce procédé de réduction de bruit pour une séquence d'images fluoroscopiques représentant un objet ayant un mouvement de translation et de rotation ;
- la figure 5 représente un exemple de sous-étapes du procédé de réduction de bruit permettant de réaliser un filtrage spatial ;
- la figure 6 représente un exemple de fonction de seuillage utilisée pour le filtrage spatial de la figure 5.
Pour la suite de la description, on considère une séquence d'images fluoroscopiques xn acquises par un détecteur à rayons X, où n désigne les différents instants d'acquisition dans la séquence. Chaque image est composée d'une matrice de pixels ayant chacun une valeur xn(m, p) représentative d'un niveau de signal, avec m et p désignant respectivement la ligne et la colonne de la matrice. La valeur est par exemple une valeur numérique comprise entre 0 et 255. La valeur numérique est ainsi codée sur 8 bits. Elle peut également être codée sur 14 bits. Le rapport contraste sur bruit (CNR) est défini comme suit :
\contraste— bruit\
CNR =
σ
où le bruit est de préférence déterminé sur une zone de l'image où il n'y a pas d'objets d'intérêt, où le contraste est de préférence déterminé à partir de points correspondant à un objet d'intérêt dans l'image, et où σ désigne l'écart type du bruit, c'est-à-dire l'écart type des valeurs de pixels correspondant à une zone de l'image ne comportant pas d'objets d'intérêt. La figure 1 représente les points considérés (en blanc) pour déterminer le contraste d'une endoprothèse vasculaire, mieux connue sous la dénomination anglo- saxonne de "stent".
La figure 2 représente des étapes possibles pour le procédé de réduction de bruit dans la séquence d'images fluoroscopiques. Le procédé est décrit pour l'image xn acquise à l'instant n. Il peut cependant être appliqué successivement sur chacune des images de la séquence.
Dans une première étape 21 , l'image xn est traitée par une chaîne de prétraitement. Le prétraitement peut notamment comporter des regroupements de pixels adjacents dans l'image. Les pixels peuvent être regroupés par exemple en fonction de la résolution souhaitée ou du niveau de bruit souhaité. Le prétraitement peut également comporter une phase d'étalonnage du détecteur à rayons X. En particulier, le gain et l'offset à appliquer à chaque pixel du détecteur peuvent être déterminés. En outre, l'étalonnage peut également comprendre une identification des pixels prenant une valeur aléatoire ou une valeur constante quel que soit le rayonnement reçu. Ces pixels défectueux peuvent être corrigés en interpolant leur valeur par les valeurs des pixels voisins. Le prétraitement peut aussi comporter une remise à l'échelle des images, c'est-à-dire un ajustement global de la valeur de tous les pixels de chaque image, de manière à ce que la valeur moyenne des pixels soit constante pour toutes les images de la séquence. Enfin, le prétraitement peut comporter une sous- étape d'adaptation des images. Étant donnée la physique des rayons X, le bruit présent dans chaque image xn est un bruit poissonien. Or un tel bruit est inadapté aux traitements de réduction de bruit. De plus, chaque image xn est une image représentative des atténuations des rayons X, alors qu'une image des épaisseurs radiologiques est plus naturelle à analyser. La formule d'Ascombe permet d'obtenir une image des épaisseurs radiologiques avec un bruit pseudo-gaussien :
Figure imgf000010_0001
où x0 n(m, p) est la valeur du pixel (m, p) avant la sous-étape d'adaptation de l'image xn.
Dans une deuxième étape 22, l'image xn est filtrée par un filtre temporel. Le filtrage temporel est par exemple un filtrage appelé "Fluoroscopic Noise Réduction" (FNR). Ce filtrage est décrit notamment dans le brevet US 6,314,160 B1 . Soit yn l'image xn obtenue après le filtrage temporel. On note A(m, p), la différence entre la valeur xn(m, p) dans l'image xn acquise à l'instant n et la valeur yn_i(m, p) dans l'image filtrée n_! correspondant à l'acquisition à l'instant n-1 :
A(m, p) = xn(m, p) - y„_i(m, p)
Pour alléger les notations, on notera simplement xn la valeur d'un pixel de l'image xn, yn la valeur d'un pixel de l'image yn et Δ la différence entre les deux valeurs xn et yn_1 du pixel considéré. Le filtrage FNR se décline selon deux modes, appelés mode adaptatif et mode non adaptatif. Le mode adaptatif est appelé ainsi du fait de sa prise en compte du mouvement des objets. Le filtrage FNR consiste à comparer pixel à pixel, la différence Δ à des seuils prédéterminés. Le premier seuil Si correspond à une différence Δ relativement faible. Le deuxième seuil S2 correspond à une différence Δ relativement forte. Le filtrage FNR est alors paramétré par des coefficients a et β en fonction du résultat de cette différence Δ et du mode considéré. On note a la force du filtre :
Figure imgf000011_0001
où σχ est l'écart type dans l'image xn et où oy est l'écart type dans l'image
On définit aussi les constantes AA, BB, CC et DD, comme suit :
. . l-a
■ A4 =
2.σχ
■ BB = 2a - 1
. cc = - a^ 2
■ DD = σ(1 - a)
Pour le mode adaptatif, les seuils Si et S2 sont définis comme suit :
' S1 = σ
' S2 = 2σ
et les coefficients a et β sont déterminés comme suit :
. , , ra = a
■ si \A\ < Sl t \β = 0
! _ AAA2 +BB.\K\+CC
A ~ Û
/? = 0
[ a = 1——
■ sinon, ΙΔΙ
1 3 = 0
Pour le mode non adaptatif, les seuils Si et S2 sont définis comme suit :
Figure imgf000011_0002
et les coefficients a et β sont déterminés comme suit :
■
Figure imgf000011_0003
(a = 1
■ sinon, \R _ « _
[p 3 V
Le filtrage FNR proprement dit consiste à déterminer une nouvelle valeur yn pour chaque pixel de l'image à partir de sa valeur xn et des coefficients a et β :
yn = a. xn + (1 - a) n_! + β
Autrement dit, la valeur xn de chaque pixel est corrigée à partir des coefficients a et β. Sur la figure 2, l'étape 22 de filtrage FNR est représentée pour le mode adaptatif. A titre d'exemple, dans une première sous-étape 221 , la valeur absolue de la différence Δ est comparée pour chaque pixel xn au seuil Sl t à savoir l'écart type o. Si la valeur absolue de la différence Δ est inférieure au seuil Sl t le pixel considéré représente probablement un niveau de bruit et un filtrage dit fort est appliqué dans une sous-étape 222 avec les paramètres correspondants. Sinon, dans une sous-étape 223, la valeur absolue de la différence Δ est comparée pour chaque pixel xn au seuil S2, à savoir 2o. Si la valeur absolue de la différence Δ est supérieure au seuil S2 , le pixel considéré correspondant probablement à un point de l'objet ayant bougé et un filtrage dit faible est appliqué dans une sous-étape 224 avec les paramètres correspondants. Enfin, si la valeur absolue de la différence Δ est comprise entre les seuils 5X et S2, une incertitude existe. Afin d'obtenir une fonction continue entre le filtrage fort et le filtrage faible, un filtrage moyen est appliqué dans une sous-étape 225 avec les paramètres correspondants.
Dans une troisième étape 23, l'image yn obtenue après filtrage temporel est filtrée par un filtre spatial basé sur la transformée en curvelets. La théorie des curvelets utilise le principe des espaces multirésolution ou multi-échelle, une échelle correspondant à un niveau de partition de l'espace. L'espace est partitionné en curvelets dont les enveloppes sont de dimensions 2j x 2j/2 , où j est un entier positif désignant l'échelle. Chaque curvelet dépend du nombre d'échelles considérées, de sa position dans l'image et de son orientation. Pour chaque échelle, la meilleure approximation possible de l'image est calculée. D'une échelle à l'autre, seules les améliorations obtenues sont représentées. Le nombre d'échelles détermine la qualité de la reconstruction de l'image. Avec les curvelets, un nombre d'échelles compris entre 6 et 10 permet de fournir des résultats très satisfaisants en termes de rapport contraste sur bruit. Dans une première sous-étape 231 de l'étape de filtrage spatial 23, l'image yn est transformée du domaine spatial au domaine des curvelets par une transformation en curvelets. De préférence, une transformation en curvelets discrète est utilisée. Une telle transformation est notamment décrite dans E. Candès, Demanet L, Donoho D. and Ying L, "Fast Discrète Curvelet Transforme", Multiscale Model. Simul., vol. 5, no. 3, 2006, pp. 861 - 899. Sous leur forme discrète, les curvelets ne peuvent pas être échantillonnées selon les grilles de l'espace de Fourier. Le document de Donoho et al propose ainsi deux solutions pour translater les curvelets à échelle et orientation données sur la grille spatiale. Une première solution pour transformer l'image dans le domaine des curvelets discrètes est appelée transformée en curvelets discrète via USFFT, où USFFT est le sigle anglo-saxon pour "unequally-spaced fast Fourier transform" (USFFT), que l'on peut traduire en français par "transformée de Fourier rapide pour données à répartition irrégulière". Dans cette solution, la grille de translation est inclinée de manière à être alignée dans la direction des curvelets à échelle et angle donnés. Une deuxième solution est appelée transformée en curvelets discrète via wrapping, où le terme "wrapping" peut être traduit en français par "enroulement". Dans cette solution, la grille de translation est identique pour tous les angles d'un même quadrant à une échelle donnée. Cela implique que le nombre d'angles utilisés pour la transformée soit un multiple de 4. La transformée en curvelets via USFFT est plus fidèle à la théorie continue des curvelets que la transformée en curvelets via wrapping. Cependant, son coût de calcul est plus important. Dans les deux solutions, la transformée en curvelets de l'image yn peut être représentée par un ensemble de coefficients c (j, l, k), où l'exposant D fait référence à la transformée discrète, et où j, l et k désignent respectivement l'échelle, l'orientation et la position de la curvelet considérée. Selon une forme préférentielle de réalisation du filtre spatial, le nombre d'échelles utilisées est égal à 9. Selon une autre forme de réalisation, le nombre d'échelles utilisées est égal à 6. Cette dernière forme de réalisation implique un nombre plus faible de coefficients à calculer et donc une durée de traitement inférieure. Par ailleurs, le nombre d'orientations utilisées est avantageusement égal à 16. Ce nombre présente un bon compromis entre la résolution en orientation et le coût de calcul des curvelets, notamment pour des images réelles dont les courbes nécessitent un grand nombre d'orientations pour les approcher. De plus, le nombre 16 étant un multiple de 4, il facilite la transformée en curvelets par wrapping.
Dans une deuxième sous-étape 232 de l'étape de filtrage spatial 23, les coefficients c (j, l, k) de l'image yn sont seuillés par une fonction de seuillage. L'étape de seuillage 232 est basée sur l'hypothèse selon laquelle certains coefficients c (j, l, k) représentent le bruit alors que d'autres contiennent l'information utile. L'étape de seuillage 232 consiste donc à identifier et à supprimer les coefficients ne contenant pas d'informations utiles. La fonction de seuillage se caractérise d'une part par sa forme et d'autre part par son seuil T. Il existe principalement deux fonctions de seuillage, la fonction de seuillage dit doux et la fonction de seuillage dit dur. En simplifiant la notation des coefficients par cn et en appelant cn les coefficients seuillés, la fonction de seuillage dur peut être modélisée comme suit :
Figure imgf000014_0001
La fonction de seuillage doux peut être modélisée comme suit :
Figure imgf000014_0002
sinon
La figure 3A représente par une courbe 31 la fonction de seuillage dur et la figure 3B représente par une courbe 32 la fonction de seuillage doux. La fonction de seuillage doux a une capacité de débruitage plus importante que la fonction de seuillage dur. En revanche, elle modifie la valeur des pixels en fonction du seuil T choisi. Il en résulte une perte de précision sur certains détails de l'image. La détermination du seuil T est une étape cruciale car elle détermine la performance du filtrage spatial, tant en ce qui concerne la capacité de débruitage que l'introduction d'artéfacts.
Selon une première forme de réalisation, le seuil T est déterminé comme étant le seuil universel Τυ décrit notamment dans D.L. Donoho and Johnstone, I.M., "Idéal spatial adaptation via wavelet shrinkage", Biometrica, vol. 81 , 1994, pp. 425-455. Ce seuil Τυ est défini par la relation suivante :
Figure imgf000014_0003
où N est le nombre d'échantillons, c'est-à-dire ici le nombre de pixels dans l'image yn. Le seuil universel Τυ n'est pas optimal et peut introduire des artéfacts. Néanmoins, des seuils proportionnels au seuil universel ry peuvent présenter de bons résultats si le facteur de proportionnalité est bien choisi.
Selon une deuxième forme de réalisation, le seuil T est déterminé par la méthode décrite dans J.L. Starck, Candès I.J. and Donoho L, "The Curvelet Transform for Image Denoising", IEEE Transactions on Image Processing, vol. 1 1 , no. 6, 2002, pp. 670-684. Cette méthode consiste à créer une image, notée Y, de mêmes dimensions que l'image yn à débruiter et de valeur uniforme 1 sur tous les pixels. L'image Y est transformée par la transformée de Fourier. La fréquence de coordonnées (0,0) dans l'image Y est centrée et son amplitude est normalisée par les dimensions de l'image. Soit F l'amplitude de la fréquence centrale. Son amplitude normalisée, notée
\F\, est par exemple réalisée selon la relation suivante :
6 F
\F\ =— ,
5 x P
où M et P désignent les dimensions de l'image Y, c'est-à-dire également de l'image xn et de l'image yn.
On obtient ainsi une image avec un pic de Dirac centré sur l'image, dont l'amplitude est fonction des dimensions de l'image. Le seuil de cette méthode, noté Ts dépend de l'échelle j et de l'orientation l. Il est déterminé en fonction de la norme des coefficients ci Y k des curvelets de l'image Y pour une échelle j et une orientation l données. Cette norme ci Y est déterminée par la relation suivante :
Figure imgf000015_0001
La fonction de seuillage est une fonction de seuillage dur, définie comme suit :
_ fO si |cn | < Ts = r. ¾. a
cn sinon
Le facteur multiplicatif r est déterminé en fonction de l'échelle j. Pour l'échelle la plus fine (par exemple j = 9), le facteur r est égal à 4. Pour les autres échelles, le facteur r est égal à 3. Ces facteurs présentent une bonne capacité de débruitage pour des images présentant un fort rapport contraste sur bruit. Cependant, pour des images fluoroscopiques, des artéfacts sont introduits. Ces artéfacts peuvent être évités en adaptant le facteur r. En particulier, pour l'échelle la plus grossière (y = 1), le facteur r peut être égal à 0. Il n'y a donc pas de seuillage. Pour l'échelle la plus fine, le facteur r peut encore être égal à 4. Et pour les autres échelles, le facteur r peut être égal à 1.
Dans une troisième sous-étape 233 de l'étape de filtrage spatial 23, l'image yn est transformée du domaine des curvelets au domaine spatial (matriciel) par une transformation en curvelets inverse. A cet effet, les coefficients seuillés sont utilisés. L'image obtenue après filtrage temporel et filtrage spatial est notée zn.
La figure 4 illustre par un graphique la performance du procédé de réduction de bruit selon l'invention pour une séquence d'images représentant un stent ayant un mouvement de translation et de rotation. Sur le graphique, l'axe des abscisses représente le CNR moyen de la séquence d'images xn et l'axe des ordonnées représente le CNR moyen de la séquence d'images obtenue après un traitement de réduction de bruit. Une première courbe 41 représente le CNR moyen de la séquence d'images zn obtenue par le procédé de réduction de bruit selon l'invention, c'est-à-dire par un filtrage temporel suivi d'un filtrage spatial par transformation en curvelets. Une deuxième courbe 42 représente le CNR moyen de la même séquence d'images après traitement par le même filtrage temporel suivi d'un filtrage spatial par transformation en ondelettes. Le graphique montre que pour une séquence d'images à traiter dont le CNR est relativement faible (0,5 - 0,7 - 1 ), le procédé de réduction de bruit selon l'invention présente une amélioration notable en termes de CNR. La réduction de bruit est encore bien meilleure pour une séquence d'images dont le CNR est plus élevé (1 ,5 - 2 - 2,5). De plus, le procédé de réduction de bruit selon l'invention permet de réduire significativement l'introduction d'artéfacts.
L'étape 23 de filtrage spatial décrite plus haut permet de réduire efficacement le niveau de bruit dans les images. Cependant, cette étape présente l'inconvénient d'introduire des artefacts apparaissant sous la forme de curvelets, c'est-à-dire des objets présentant des formes allongées et ondulantes. En effet, ces artéfacts sont dus à des coefficients ayant une grande amplitude, mais ne représentant que du bruit. Ils ne sont présents que lorsque le niveau de bruit est élevé, ce qui est le cas des images fluoroscopiques. Une solution consiste à fixer des seuils relativement bas, ce qui réduit les artefacts mais limite les performances en termes de réduction de bruit. La suite de la description décrit une autre méthode pour réaliser l'étape de filtrage spatial permettant de réduire significativement la présence d'artefacts sans réduire les performances du filtrage. Cette méthode repose sur la corrélation qui existe entre un coefficient donné et les coefficients qui lui sont liés, à savoir son coefficient parent, ses coefficients voisins et ses coefficients cousins. En particulier, il est fait l'hypothèse que si un coefficient représente du signal utile, les coefficients qui lui sont liés ont une forte probabilité de représenter également du signal utile. À l'inverse, si le module d'un coefficient est important en raison de la présence d'un artéfact, il est peu probable que les coefficients qui lui sont liés soient aussi impactés par cet artefact.
La figure 5 représente un autre exemple de sous-étapes permettant de réaliser l'étape 23 de filtrage spatial. Dans une première sous- étape 51 , l'image yn obtenue après filtrage temporel est transformée du domaine spatial au domaine des curvelets. La transformée en curvelets de l'image yn est également représentée par des coefficients notés c° (j,l,k) , ou cn sous forme simplifiée. De préférence, la transformation en curvelets est faite avec 6 échelles et 12 orientations. Selon une forme particulière de réalisation, l'image yn est représentée par des ondelettes pour l'échelle la plus fine, par exemple pour l'échelle j = 5 .
Dans une deuxième sous-étape 52, trois seuils notés Tx , T2 et T3 sont déterminés tels que Tl < T2 < T3 . Ces seuils sont déterminés en fonction de l'écart-type asfl de la sous-bande des coefficients. La sous-bande d'un coefficient donné c° (j,l,k) regroupe l'ensemble des coefficients ayant la même échelle j et la même orientation / . Il s'agit donc de l'ensemble des positions pour l'échelle et l'orientation considérées. Les seuils Τλ , T2 et T3 définissent une première zone pour les valeurs inférieures au seuil Τλ , une deuxième zone pour les valeurs comprises entre les seuils Tx et T2 , une troisième zone pour les valeurs comprises entre les seuils T2 et T3 , et une quatrième zone pour les valeurs supérieures au seuil T3 . À titre d'exemple, on considère les valeurs de seuil suivantes : Tx = 7,5 , T2 = \ 1,25 , et T3 = 15 .
Dans une troisième sous-étape 53, les coefficients sont corrigés en fonction de leur coefficient parent et de leurs coefficients voisins. De préférence, on considère uniquement les coefficients voisins non corrélés. On note c°(j,l,k) , ou cn sous forme simplifiée, les coefficients corrigés.
Pour un coefficient donné c°(j,l,k), son coefficient parent est le coefficient correspondant à l'échelle plus grossière précédente, c'est-à-dire le coefficient c°(j-l,l,k). Ses coefficients voisins sont les coefficients c°(j,l,k') correspondant aux positions voisines k' , c'est-à-dire les positions contiguës ayant la même échelle et la même orientation. Pour une position repérée par le couple {x,y), il s'agit donc des positions (x + \,y), (x-\,y), (x,y + \), (x,y-l), (JC + 1,J> + I), (JC + 1,J>-1),(JC-1,J> + I) et (jc-l,j>-l). On qualifie de "corrélés" les coefficients voisins situés dans la direction de la sous-bande. Pour chaque coefficient c°(j,l,k), le module de son coefficient parent et la moyenne des modules de ses coefficients voisins non corrélés sont comparés aux seuils Τλ , T2 et T3. Trois cas sont à considérer. Le premier cas est celui où le module du coefficient considéré c°(j,l,k) se situe soit dans la même zone que le module de son coefficient parent c°(j-l,l,k), soit dans la même zone que la moyenne des modules de ses coefficients voisins c°(j,l,k'). Dans ce premier cas, le module du coefficient considéré reste inchangé (c°(j,l,k) = c°(j,l,k)). Le deuxième cas est celui où le module du coefficient parent c°(j-l,l,k) et la moyenne des modules des coefficients voisins c°(j,l,k') se situent dans la même zone, différente de celle du module du coefficient considéré c (/,/,&). Le module du coefficient corrigé c°(j,l,k) prend alors la valeur minimale de la zone dans laquelle se situent le module du coefficient parent et la moyenne des modules des coefficients voisins si le module du coefficient considéré c°(j,l,k) se situe dans une zone à plus faibles valeurs, et la valeur maximale de cette zone sinon. Le troisième cas est celui où le module du coefficient considéré c°(j,l,k), le module de son coefficient parent c°(j-l,l,k), et la moyenne des modules de ses coefficients voisins c°(j,l,k') se situent dans trois zones distinctes. Le module du coefficient corrigé c (j,l,k) prend alors la valeur minimale de la zone dans laquelle se situe la moyenne des modules des coefficients voisins si le module du coefficient considéré c°(j,l,k) se situe dans une zone à plus faibles valeurs, et la valeur maximale de cette zone sinon. Le troisième cas privilégie la zone des coefficients voisins non corrélés au détriment de la zone du coefficient parent du fait qu'à l'échelle la plus grossière, les coefficients ne possèdent pas de parent. Il serait néanmoins envisageable de corriger les coefficients des autres échelles en fonction de la zone dans laquelle se situe le module de leur coefficient parent.
Dans une quatrième sous-étape 54, les coefficients corrigés cn sont seuillés par une fonction de seuillage. On note c° (j,l,k) , ou cn sous forme simplifiée, les coefficients corrigés seuillés. La figure 6 représente la fonction de seuillage correspondant aux valeurs de seuils Tx , T2 et T3 considérées. Cette fonction est modélisée comme suit :
Figure imgf000019_0001
où f{cn ) est une sigmoïde définie comme suit :
f(cn ) = L_
1 + exp(- A(c„ - m))
Les paramètres m et À sont déterminés de manière à ce que la fonction de seuillage soit continue. Ils sont ainsi définis comme suit :
Figure imgf000019_0002
- ln(( ).0l)
T - - m
La fonction de seuillage annule donc les coefficients dont le module est compris dans la première zone 61 , elle applique un seuillage doux aux coefficients dont le module est compris dans la deuxième zone 62, elle applique un seuillage dur aux coefficients dont le module est compris dans la quatrième zone 64, et elle applique un seuillage de transition aux coefficients dont le module est compris dans la troisième zone 63. Le seuillage doux permet de réduire le niveau de bruit, mais il étale les contours. Le seuillage dur conserve bien les contours mais est moins efficace en termes de réduction de bruit. L'utilisation conjointe de ces deux types de seuillage permet de combiner leurs avantages respectifs en fonction des zones. Lorsque l'image yn est représentée par des ondelettes pour l'échelle la plus fine, les coefficients cn sont seuillés par une fonction de seuillage dur modélisée comme suit :
Figure imgf000020_0001
où le seuil T4 est défini à partir de l'écart-type asfl de la sous-bande du coefficient c„ considéré :
4 χ
Dans une cinquième sous-étape 55, l'image yn est transformée du domaine des curvelets au domaine spatial par une transformation en curvelets inverse. À cet effet, les coefficients seuillés corrigés sont utilisés. L'image obtenue après filtrage temporel et filtrage spatial est notée z'n .
Le filtrage spatial décrit en référence à la figure 5 utilise à la fois les coefficients parents et les coefficients voisins non corrélés. Le filtrage spatial pourrait également utiliser les coefficients voisins corrélés, ou l'ensemble des coefficients voisins. Il pourrait aussi utiliser uniquement les coefficients parents ou les coefficients voisins, ou encore combiner ces coefficients avec les coefficients cousins. Ces coefficients sont ceux ayant la même échelle et la même position que le coefficient de référence. Le filtrage spatial pourrait également ne pas comporter de sous-étape de correction des coefficients, la sous-étape 54 de seuillage étant dans ce cas réalisée directement sur les coefficients cn des transformées en curvelets.

Claims

REVENDICATIONS
1 . Procédé de réduction du bruit dans une séquence d'images fluoroscopiques acquises par un détecteur à rayons X, chaque image (xn) étant formée dans le domaine spatial d'une matrice de pixels, chaque pixel ayant une valeur représentative d'un niveau de signal, le procédé comprenant les étapes successives suivantes pour chaque image (xn) :
■ appliquer un filtrage temporel (22) sur l'image (xn) acquise à un instant n, l'application dudit filtrage temporel comportant les sous-étapes suivantes pour chaque pixel de l'image :
- déterminer (221 , 223) une différence entre la valeur du pixel considéré dans l'image (xn) acquise à l'instant n et la valeur du pixel correspondant dans l'image (yn_i) acquise à un instant n-1 sur laquelle le filtrage temporel (22) a été appliqué,
- si la différence est inférieure à un premier seuil Si prédéterminé, corriger (222) la valeur du pixel correspondant avec un premier ensemble de paramètres,
- si la différence est supérieure à un deuxième seuil S2 prédéterminé, corriger (224) la valeur du pixel correspondant avec un deuxième ensemble de paramètres,
■ appliquer un filtrage spatial (23) sur l'image (yn) acquise à l'instant n, l'application dudit filtrage spatial comportant les sous-étapes suivantes :
- transformer (231 , 51 ) l'image (yn) acquise à l'instant n du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image dans le domaine des curvelets étant représentée par un ensemble de coefficients,
- seuiller (232, 54) les coefficients de l'image (yn) par une fonction de seuillage, la fonction de seuillage annulant les coefficients inférieurs à un troisième seuil prédéterminé, et conservant ou ajustant les coefficients supérieurs au troisième seuil prédéterminé,
- transformer (233, 55) l'image (yn) dont les coefficients ont été seuillés du domaine des curvelets au domaine spatial par une transformation en curvelets inverse.
2. Procédé selon la revendication 1 , dans lequel l'étape (231 ) de transformation de l'image (yn) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets à 6 échelles.
3. Procédé selon la revendication 1 , dans lequel l'étape (231 ) de transformation de l'image (yn) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets à 9 échelles.
4. Procédé selon l'une des revendications précédentes, dans lequel l'étape (231 ) de transformation de l'image (yn) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets comprenant 16 orientations.
5. Procédé selon l'une des revendications précédentes, dans lequel l'étape (231 ) de transformation de l'image (yn) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets discrète par USFFT.
6. Procédé selon l'une des revendications 1 à 4, dans lequel l'étape (231 ) de transformation de l'image (yn) acquise à l'instant n dans le domaine des curvelets utilise une transformée en curvelets discrète par wrapping.
7. Procédé selon l'une des revendications précédentes, dans lequel la fonction de seuillage est une fonction de seuillage dur, c'est-à-dire que les coefficients sont annulés s'ils sont inférieurs au troisième seuil prédéterminé, et conservés sinon.
8. Procédé selon la revendication 7, dans lequel le troisième seuil prédéterminé Ts est déterminé par les étapes successives suivantes :
■ créer une image Y de mêmes dimensions que l'image (yn) acquise à l'instant n et de valeur uniforme 1 sur tous les pixels,
■ transformer l'image Y par une transformée de Fourier,
■ normaliser l'amplitude F du pic dans la transformée de Fourier de l'image Y par la relation : F
5 x P
où M et P désignent les dimensions de l'image Y,
■ transformer la transformée de Fourier de l'image Y du domaine spatial au domaine des curvelets par une transformation en curvelets, chaque image dans le domaine des curvelets étant représentée par un ensemble de coefficients cj k,
■ déterminer une norme cjtl des coefficients cjl k par la relation :
Figure imgf000023_0001
où M et P désignent les dimensions de l'image Y et / désigne les orientations des curvelets,
■ déterminer le seuil prédéterminé Ts par la relation :
Ts = r. cj o
où σ est l'écart type du bruit dans l'image (yn) acquise à l'instant n et sur laquelle le filtrage temporel (22) a été appliqué, et où le facteur r prend la valeur 0 pour l'échelle la plus grossière des curvelets, la valeur 4 pour l'échelle la plus fine, et la valeur 1 pour les autres échelles.
9. Procédé selon l'une des revendications 1 à 6, dans lequel la fonction de seuillage applique un seuillage doux aux coefficients dont le module est compris entre le troisième seuil prédéterminé Ti et un quatrième seuil prédéterminé T2, un seuillage dur aux coefficients dont le module est supérieur à un cinquième seuil prédéterminé T3, et un seuillage de transition entre le seuillage doux et le seuillage dur aux coefficients dont le module est compris entre le quatrième seuil prédéterminé T2 et le cinquième seuil prédéterminé T3, les seuils prédéterminés T1 , T2, et T3 étant tels que T1 < T2 < T3.
10. Procédé selon la revendication 9, dans lequel le seuillage de transition est définie par une sigmoïde définie comme suit :
Figure imgf000023_0002
où cn sont les coefficients devant être seuillés, et où les paramètres m et À sont déterminés comme suit :
Figure imgf000024_0001
11 . Procédé selon l'une des revendications 9 et 10, dans lequel l'étape (23) d'application du filtrage spatial comprend, en outre, préalablement à la sous-étape (54) de seuillage :
■ une sous-étape (53) de correction des coefficients en fonction de leur coefficient parent, de leurs coefficients voisins, et/ou de leurs coefficients cousins, le coefficient parent d'un coefficient donné c° (j, l, k) étant le coefficient c° (j - l, l, k) ayant la même orientation et la même position que le coefficient donné à l'échelle plus grossière précédente j - 1, ses coefficients voisins étant les coefficients c° (j, l, k') correspondant aux positions contiguës de la position du coefficient donné et ayant la même échelle et la même orientation, les coefficients cousins étant les coefficients ayant la même échelle et la même position que le coefficient donné.
12. Procédé selon la revendication 11 , dans lequel une première zone de valeurs (61 ) est définie pour les valeurs inférieures au troisième seuil prédéterminé ΤΊ, une deuxième zone de valeurs (62) est définie entre le troisième seuil prédéterminé ΤΊ et le quatrième seuil prédéterminé T2, une troisième zone de valeurs (63) est définie entre le quatrième seuil prédéterminé T2 et le cinquième seuil prédéterminé T3, et une quatrième zone de valeurs (64) est définie pour les valeurs supérieures au cinquième seuil prédéterminé T3, la sous-étape (54) de correction des coefficients appliquant les corrections suivantes :
■ lorsque le module du coefficient considéré c° (j, l, k) se situe soit dans la même zone que le module de son coefficient parent c° (j - l, l, k) , soit dans la même zone que la moyenne des modules de ses coefficients voisins c (/,/,&'), le module du coefficient corrigé c (j, l, k) prend la valeur du module du coefficient c° (j, l, k) , ■ lorsque le module du coefficient parent c° {j - l, l, k) et la moyenne des modules des coefficients voisins c° (j, l, k') se situent dans une même zone, différente de celle du module du coefficient considéré c° (j, l, k) , le module du coefficient corrigé c (j, l, k) prend la valeur minimale de la zone dans laquelle se situent le module du coefficient parent et la moyenne des modules des coefficients voisins si le module du coefficient considéré c° (j, l, k) se situe dans une zone à plus faibles valeurs, et la valeur maximale de cette zone sinon,
■ lorsque le module du coefficient considéré c (/,/,&), le module de son coefficient parent c° (j - l, l, k) , et la moyenne des modules de ses coefficients voisins c° (j, l, k') se situent dans trois zones distinctes, le module du coefficient corrigé c {j, l, k) prend la valeur minimale de la zone dans laquelle se situe la moyenne des modules des coefficients voisins si le module du coefficient considéré c (/,/,&) se situe dans une zone à plus faibles valeurs, et la valeur maximale de cette zone sinon.
13. Procédé selon l'une des revendications précédentes, dans lequel l'étape (22) d'application du filtrage temporel sur l'image acquise (xn) à l'instant n comprend, à la suite de la sous-étape (221 , 223) de détermination de la différence entre la valeur du pixel considéré dans l'image acquise (xn) à l'instant n et la valeur du pixel correspondant dans l'image (yn_i) acquise à l'instant n-1 sur laquelle le filtrage temporel (22) a été appliqué, la sous-étape supplémentaire suivante :
■ si la différence est comprise entre le premier seuil Si prédéterminé et le deuxième seuil S2 prédéterminé, corriger (225) la valeur du pixel correspondant avec un troisième ensemble de paramètres.
PCT/EP2012/064347 2011-07-22 2012-07-20 Procédé de réduction du bruit dans une séquence d'images fluoroscopiques par filtrage temporel et spatial WO2013014110A1 (fr)

Priority Applications (3)

Application Number Priority Date Filing Date Title
US14/234,292 US9058654B2 (en) 2011-07-22 2012-07-20 Method for reducing noise in a sequence of fluoroscopic images by temporal and spatial filtering
EP12737811.5A EP2734975B1 (fr) 2011-07-22 2012-07-20 Procédé de réduction du bruit dans une séquence d'images fluoroscopiques par filtrage temporel et spatial
JP2014520683A JP6120843B2 (ja) 2011-07-22 2012-07-20 時間フィルタリングおよび空間フィルタリングによる一連の蛍光透視画像における雑音を低減させるための方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
FR1102306A FR2978273B1 (fr) 2011-07-22 2011-07-22 Procede de reduction du bruit dans une sequence d'images fluoroscopiques par filtrage temporel et spatial
FR11/02306 2011-07-22

Publications (1)

Publication Number Publication Date
WO2013014110A1 true WO2013014110A1 (fr) 2013-01-31

Family

ID=46548470

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/EP2012/064347 WO2013014110A1 (fr) 2011-07-22 2012-07-20 Procédé de réduction du bruit dans une séquence d'images fluoroscopiques par filtrage temporel et spatial

Country Status (5)

Country Link
US (1) US9058654B2 (fr)
EP (1) EP2734975B1 (fr)
JP (1) JP6120843B2 (fr)
FR (1) FR2978273B1 (fr)
WO (1) WO2013014110A1 (fr)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2016527935A (ja) * 2013-06-28 2016-09-15 タレス 蛍光透視画像シーケンスにおけるノイズを減らす方法

Families Citing this family (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2654015A1 (fr) * 2012-04-21 2013-10-23 General Electric Company Procédé, système et support lisible sur ordinateur pour le traitement dýune image vidéo médicale
JP6381198B2 (ja) * 2013-11-08 2018-08-29 キヤノン株式会社 制御装置、制御方法及びプログラム
FR3018147B1 (fr) * 2014-03-03 2016-03-04 Sagem Defense Securite Debruitage video optimise pour systeme multicapteur heterogene
CN104463325A (zh) * 2014-12-04 2015-03-25 中国科学院电子学研究所 一种极地探冰雷达原始数据噪声抑制方法
CN105259539B (zh) * 2015-10-29 2017-09-29 中国科学院大学 一种针对极地探冰雷达数据的分部式噪声抑制方法
US9934557B2 (en) 2016-03-22 2018-04-03 Samsung Electronics Co., Ltd Method and apparatus of image representation and processing for dynamic vision sensor
US11113791B2 (en) * 2017-01-03 2021-09-07 Flir Systems, Inc. Image noise reduction using spectral transforms
US10983236B2 (en) 2017-06-20 2021-04-20 Saudi Arabian Oil Company Super-resolution radon transform based on thresholding
JP7183746B2 (ja) * 2018-12-03 2022-12-06 コニカミノルタ株式会社 制御装置及びプログラム
JP7257829B2 (ja) * 2019-03-15 2023-04-14 ソニー・オリンパスメディカルソリューションズ株式会社 画像処理装置、画像処理方法およびプログラム
CN112051612A (zh) * 2019-06-05 2020-12-08 河北地质大学 一种消除地震数据随机噪声的方法和系统
DE102019122667A1 (de) * 2019-08-22 2021-02-25 Schölly Fiberoptic GmbH Verfahren zur Unterdrückung von Bildrauschen in einem Videobildstrom, sowie zugehöriges medizinisches Bildaufnahmesystem und Computerprogrammprodukt
CN111766306B (zh) * 2020-05-21 2023-06-27 广州多浦乐电子科技股份有限公司 一种全聚焦成像方法及装置
CN112633202B (zh) * 2020-12-29 2022-09-16 河南大学 一种基于双重去噪联合多尺度超像素降维的高光谱图像分类算法
CN112800727B (zh) * 2021-04-14 2021-07-20 北京三维天地科技股份有限公司 给pdf文件加批注的方法及应用系统
CN116883292B (zh) * 2023-09-07 2023-11-28 上海海栎创科技股份有限公司 基于伪高斯的图像方向场获取方法、系统及计算机设备

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6314160B1 (en) 1999-12-17 2001-11-06 General Electric Company Method and apparatus for performing fluoroscopic noise reduction
US7542622B1 (en) * 2003-06-02 2009-06-02 The Trustees Of Columbia University In The City Of New York Spatio-temporal treatment of noisy images using brushlets
WO2009100032A1 (fr) * 2008-02-05 2009-08-13 Ntt Docomo, Inc. Diminution du bruit et/ou du scintillement dans des séquences vidéo au moyen du traitement spatial et temporel
EP2151193A1 (fr) * 2007-05-09 2010-02-10 Hitachi Medical Corporation Dispositif à ultrasons et procédé d'amélioration de la qualité d'image du dispositif à ultrasons

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH06169906A (ja) * 1992-12-10 1994-06-21 Toshiba Corp Ｘ線画像診断装置
JP3326885B2 (ja) * 1993-06-30 2002-09-24 株式会社島津製作所 Ｘ線テレビジョン装置
DE69824230T2 (de) * 1997-04-08 2005-07-07 Koninklijke Philips Electronics N.V. Verarbeitungssystem einer verrauschten Bildsequenz und medizinisches Untersuchungsgerät mit einem solchen System
US6335990B1 (en) * 1997-07-03 2002-01-01 Cisco Technology, Inc. System and method for spatial temporal-filtering for improving compressed digital video
US6281942B1 (en) * 1997-08-11 2001-08-28 Microsoft Corporation Spatial and temporal filtering mechanism for digital motion video signals
US6018600A (en) * 1998-04-15 2000-01-25 Arch Development Corp. Sampling and reconstruction of signals and images including MR images of multiple regions
US6907143B2 (en) * 2001-05-16 2005-06-14 Tektronix, Inc. Adaptive spatio-temporal filter for human vision system models
US7840625B2 (en) * 2005-04-07 2010-11-23 California Institute Of Technology Methods for performing fast discrete curvelet transforms of data
US20100061656A1 (en) * 2006-11-08 2010-03-11 Koninklijke Philips Electronics N. V. Noise reduction of an image signal
US8553783B2 (en) * 2008-12-29 2013-10-08 Stmicroelectronics Asia Pacific Pte Ltd. Apparatus and method of motion detection for temporal mosquito noise reduction in video sequences
US8942454B2 (en) * 2009-12-03 2015-01-27 The United States Of America, As Represented By The Secretary, Department Of Health And Human Services Signal to-noise enhancement in imaging applications using a time-series of images

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6314160B1 (en) 1999-12-17 2001-11-06 General Electric Company Method and apparatus for performing fluoroscopic noise reduction
US7542622B1 (en) * 2003-06-02 2009-06-02 The Trustees Of Columbia University In The City Of New York Spatio-temporal treatment of noisy images using brushlets
EP2151193A1 (fr) * 2007-05-09 2010-02-10 Hitachi Medical Corporation Dispositif à ultrasons et procédé d'amélioration de la qualité d'image du dispositif à ultrasons
WO2009100032A1 (fr) * 2008-02-05 2009-08-13 Ntt Docomo, Inc. Diminution du bruit et/ou du scintillement dans des séquences vidéo au moyen du traitement spatial et temporel

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
D.L. DONOHO; JOHNSTONE, I.M.: "Ideal spatial adaptation via wavelet shrinkage", BIOMETRICA, vol. 81, 1994, pages 425 - 455
E. CANDÈS; DEMANET L.; DONOHO D.; YING L.: "Fast Discrete Curvelet Transforms", MULTISCALE MODEL. SIMUL., vol. 5, no. 3, 2006, pages 861 - 899
J.L. STARCK; CANDÈS I.J.; DONOHO L.: "The Curvelet Transform for Image Denoising", IEEE TRANSACTIONS ON IMAGE PROCESSING, vol. 11, no. 6, 2002, pages 670 - 684, XP002485325, DOI: doi:10.1109/TIP.2002.1014998
M. TOMIC ET AL: "Adaptive spatio-temporal denoising of fluoroscopic X-ray sequences", BIOMEDICAL SIGNAL PROCESSING AND CONTROL, vol. 7, no. 2, 5 March 2011 (2011-03-05), pages 173 - 179, XP055021074, ISSN: 1746-8094, DOI: 10.1016/j.bspc.2011.02.003 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2016527935A (ja) * 2013-06-28 2016-09-15 タレス 蛍光透視画像シーケンスにおけるノイズを減らす方法

Also Published As

Publication number Publication date
JP6120843B2 (ja) 2017-04-26
JP2014520641A (ja) 2014-08-25
US20140153842A1 (en) 2014-06-05
EP2734975B1 (fr) 2015-06-24
FR2978273B1 (fr) 2013-08-09
FR2978273A1 (fr) 2013-01-25
US9058654B2 (en) 2015-06-16
EP2734975A1 (fr) 2014-05-28

Similar Documents

Publication Publication Date Title
EP2734975B1 (fr) Procédé de réduction du bruit dans une séquence d&#39;images fluoroscopiques par filtrage temporel et spatial
CN111081354B (zh) 用于通过深度学习网络对医疗图像进行去噪的系统和方法
EP0379399B1 (fr) Procédé de calcul et d&#39;exploitation de l&#39;image en projection conique, par exemple au sens des rayons x, d&#39;un objet tridimensionnel echantillonné, et procédé de reconstruction tridimensionnelle d&#39;un objet étudié utilisant ce procédé de calcul
US7362845B2 (en) Method and apparatus of global de-noising for cone beam and fan beam CT imaging
FR2924254A1 (fr) Procede de traitement d&#39;images en radioscopie interventionnelle
JP5416912B2 (ja) データ処理装置および医用診断装置
EP3014570B1 (fr) Procede de reduction de bruit dans des sequences d&#39;images fluoroscopiques
FR2848007A1 (fr) Procede et systeme pour accentuer une image de tomosynthese en utilisant un filtrage transversal
CN102663695A (zh) 基于小波变换的dr图像去噪方法及系统
US7310436B2 (en) Systems, methods and apparatus for specialized filtered back-projection reconstruction for digital tomosynthesis
Chao et al. Astronomical image restoration using an improved anisotropic diffusion
Justen et al. A non-iterative regularization approach to blind deconvolution
Khare et al. Daubechies complex wavelet transform based multilevel shrinkage for deblurring of medical images in presence of noise
FR2848806A1 (fr) Procede de calibrage d&#39;un appareil d&#39;imagerie radiologique necessitant un nombre limite d&#39;acquisitions
Al-amri et al. A comparative study for deblured average blurred images
Lew et al. Adaptive Gaussian Wiener Filter for CT-Scan Images with Gaussian Noise Variance
Yang et al. Image deblurring using empirical Wiener filter in the curvelet domain and joint non-local means filter in the spatial domain
Saoji et al. Speckle and rician noise removal from medical images and Ultrasound images
FR2835949A1 (fr) Procede de tomosynthese a reconstruction multiplan
Bahendwar et al. A comparative performance analysis of discrete wavelet transforms for denoising of medical images
Kadri et al. Power shrinkage—curvelet domain image denoising using a new scale-dependent shrinkage function
Leung et al. Effective use of low resolution images for super-resolution reconstruction
FR2918486A1 (fr) Procede de traitement d&#39;images en radioscopie interventionnelle pour une detection de materiels d&#39;instrumentations de guidage
Bertero et al. High‐resolution imaging by multiple‐image deconvolution
CN110992436A (zh) 图像重建方法及设备、图像数据的处理方法及装置

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 12737811

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2012737811

Country of ref document: EP

ENP Entry into the national phase

Ref document number: 2014520683

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 14234292

Country of ref document: US