MULTI-QUADRATIC DYNAMIC PROGRAMMING PROCEDURE OF EDGE– PRESERVING DENOISING FOR MEDICAL IMAGES

S : In this paper, we present a computationally efficient technique for edge preserving in medical image smoothing, which is developed on the basis of dynamic programming multi-quadratic procedure. Additionally, we propose a new non-convex type of pair-wise potential functions, allow more flexibility to set a priori preferences, using different penalties for various ranges of differences between the values of adjacent image elements. The procedure of image analysis, based on the new data models, significantly expands the class of applied problems, and can take into account the presence of heterogeneities and discontinuities in the source data, while retaining high computational efficiency of the dynamic programming procedure and Kalman filter-interpolator. Comparative study shows, that our algorithm has high accuracy to speed ratio, especially in the case of high-resolution medical images.


INTRODUCTION
Low-level image processing is a functional step in almost every medical image analysis system.It is a prerequisite for proper data interpretation, diagnosis and suggestion the corresponding treatments.Medical images (example: Magnetic Resonance, Ultrasound, Computed Tomography, X-Ray), may be corrupted by a disruptive noise during acquisition and transmission process and the essential requirement for every noise reduction procedure is to preserve local image features for an accurate and effective diagnosis.Many edge-preserving denoising methods for medical images have been proposed in the literature.For example, Perona and Malik introduced nonlinear anisotropic diffusion (Perona and Malik, 1990;Gerig et al., 1992), which is an efficient method for noise reduction in MR images with local features preserving.Portilla (Portilla et al., 2003) have suggested alternative approach based on Wavelet transform.Furthermore, some other methods have good results of edge-preserving denoising: Wiener filter (Brown and Hwang, 1996), bilateral filter (Tomasi and Manduchi, 1998;Chaudhury et al., 2011), non -linear total variation (Rudin et al., 1992;Drapaca, 2009).Nevertheless, the monitoring of dynamic processes needs to improve performance of the low-level processing stage in both speed and accuracy.In this paper, we develop new non-iterative parametric procedure for edge preserving in image smoothing.This procedure can effectively remove Gaussian noise as well as Rician noise (Dekker and Sijbers, 2014), typical for MR images, with high quality.We use here Bayesian framework as one of the most popular approaches to image processing.Under this approach, the problem of image analysis can be expressed as the problem of estimation of a hidden Markov component of a two-component random field, where the analyzed image plays the role of the observed component.An equivalent representation of Markov random fields in the form of Gibbs random fields, according to Hammersley-Clifford Theorem (Hammersley and Clifford, 1971), can be used to define prior probability properties of a hidden Markov field by means of socalled Gibbs potentials for cliques.In the case of a singular loss function, the Bayesian estimation of the hidden component can be found as a maximum a posteriori probability (MAP) estimation by the minimization of the objective function often called the Gibbs energy function.The quadratic form of the Gibbs potentials corresponds to the assumption of a normal priori distribution.In the paper (Kopylov, 2010), the more general case is considered, when pair-wise Gibbs potentials are selected as a minimum of a finite set of quadratic functions, and an efficient optimization procedure is proposed for the Blake and Zisserman function (Blake and Zisserman, 1987), which is often used in edge-preserving procedures.This highly effective global optimization procedure, is based on a recurrent decomposition of the initial problem of minimizing function of many variables into a succession of partial problems, each of which consists in minimizing a function of only one variable.Such a procedure is nothing else than the enhanced version of the famous dynamic programming procedure (Kalman and Bucy, 1961).In the case of a minimum of a finite set of quadratic functions pair-wise Gibbs potentials, the procedure breaks down at each step into several parallel procedures, according to the number of quadratic functions forming the intermediate optimization problems of one variable.The corresponding intermediate objective functions are called Bellman functions and the procedure itself -a multi quadratic dynamic programming procedure.It was proven, that in the case of signal processing the number of quadratic functions that are required for representation of a Bellman function, generally, does not increase by more than one at each step (Kopylov, 2010).However, using this approach for processing the twodimensional data on the basis of tree-like approximation of the lattice neighborhood graph (Mottl, 1998), the number of quadratic functions in the Bellman functions may be too large and leads to a lack of effective implementation of the procedure.The basic idea of the proposed procedure is to find the groups of closest, in the appropriate sense, quadratic functions using k-means clustering algorithm (MacQueen, 1967), and to replace each of these groups by one quadratic function having the smallest minimun value.During the experimental study, we compare the performance of MR images denoising algorithms by using the Mean Structure Similarity Index-MSSIM (Zhou Wang, 2004 ), Signal to Noise Ratio -SNR, and Peak to Signal Noise Ratio-PSNR in the presence of Gaussian and Rician noise.Wiener filtering (WF), non -linear Total Variation (TV), Anisotropic Diffusion Filter (ADF), Fast Bilateral Filter (FBF), the Bayesian least squares with Gaussians Scale Mixture (BLS-GSM) and proposed algorithms were tested.Results show, that our algorithm has the best accuracy to speed ratio, especially in the case of highresolution images.

BAYESIAN FRAMEWORK FOR IMAGE DENOISING
Under this approach, the problem of image analysis can be expressed as the problem of estimating a hidden Markov component ( , plays the role of the observed component.An equivalent representation of Markov random fields in the form of Gibbs random fields, according to Hammersley Clifford Theorem (Hammersley JM, Clifford PE, 1971), can be used to define a priori probability properties of a hidden Markov field by means of so-called Gibbs potentials for cliques (Geman S., German D. 1984).In the case of a singular loss function, the Bayesian estimation of the hidden components can be found as Maximum a Posteriori Probability (MAP) by performing the minimization of the objective function, often called the Gibbs energy function (Mottl, 1998;Kopylov, 2010): The structure of the objective function (1) reflects the ordering property of analyzed data and determined by means of an undirected graph neighborhood G T T   .In image analysis objective function can often be represented as the sum of the two types of potential functions for cliques, called node functions and edge functions.In image smoothing node functions ( | ) x Y  t t t play the role of a penalty on the difference between the values of input data Y and the seeking function X , and are usually chosen in quadratic form imposes penalty upon the difference of values in the corresponding pair of adjacent elements of the edge ( , )   t t of neighborhood graph G , and can have various forms.The quadratic form of Gibbs potentials corresponds to the assumption of normal priori distribution Х .In a Bayesian framework, the edge functions are usually called a pair-wise potential functions.There are different pair-wise potentials in the literature.A variety of non-convex pair-wise potential functions were considered for preserving large differences of values in the corresponding pairs of adjacent elements of the edges and smoothing the other differences (Nikolova, 2010).Convex edge-preserving pair-wise potential functions were proposed to avoid the numerical involutions arising with non-convex regularization (Xiaoju and Li, 2009).However, non-convex regularization offer the best possible quality of image reconstruction with neat and exact edges (Nikolova, 2010).One of the main problems in these approaches is high computational complexity of corresponding minimization procedures which can hardly be applied to high-resolution images.In the paper (Kopylov, 2010), the different case is considered, when pair-wise Gibbs potentials are selected as a minimum of a finite set of quadratic functions, and an efficient optimization procedure for signal processing is proposed for the Blake and , which is often used in edge-preserving procedures.This highly effective global optimization procedure, is based on a recurrent decomposition of the initial problem of minimizing function of many variables into a succession of partial problems, each of which consists in minimizing a function of only one variable.Such a procedure is nothing else than one of versions of the famous dynamic programming procedure.In the case of a minimum of a finite set of quadratic functions pair-wise Gibbs potentials, the procedure breaks down at each step into several parallel procedures, according to the number of quadratic functions forming the intermediate optimization problems of one variable.The corresponding intermediate objective functions are called Bellman functions and the procedure itselfa multi quadratic dynamic programming procedure.According to the central idea of the variation approach to the analysis of ordered data sets, the result of analysis ˆ( ) X Y is defined by the condition: ( 2 ) The procedure of dynamic programming will search for optimization of the objective function (1) in the forward and backward recurrent relation.In the forward recurrent with upward recurrent relation defined by the Bellman function (Mottl, 1998): The global minimum of the Bellman function from the last variable min ( ) coincides with the global minimum of the total objective function in all variables, so the optimal value of the last variable minimize its function Bellman (Mottl, 1998): The other variables can be found by applying backward recurrent relation ( 1,.., 2,1   t N ):

MULTI-QUADRATIC DYNAMIC PROGRAMMING PROCEDURE FOR IMAGE PROCESSING
In this paper, we propose a new non-convex type of pair-wise potential functions, allows more flexibility to set a priori preferences, using different coefficients of penalty for various ranges of differences between the values of adjacent image elements: The developed image analysis procedure can significantly extend the class to solve applied problems, to take into account the presence of heterogeneities and discontinuities in the original data, while retaining the high computational efficiency of procedures of dynamic programming and Kalman filterinterpolator (Kalman and Bucy, 1961).
It can be proven, that if the pair-wise Gibbs potentials are selected as a minimum of a finite set of quadratic functions, and node functions are in quadratic form, the procedure breaks down at each step into several parallel procedures, according to the number of quadratic functions forming the intermediate optimization problems of one variable.The Bellman function at each step of the dynamic programming will has a minimum of a finite set of quadratic functions: 1 2 ( ) min ( ), ( ),..., ( ) where 2 ( ) ( ) , 0, 1,.., ( ) min ( ), ( ), , ( ) The backward recurrent relation has the following form: where , 1.. ; i j K i j   The coordinates of the intersection points on the real axis defined by the following relations: The points of intersection coordinates, if expression  under the square root is greater than or equal to zero: After that, select functions that have points of intersection with satisfying . Discard all the functions for which there is no intersection.Among the tested functions, select the function with the smallest minimum and leave it.Check its intersection with other functions.Discard all functions for which there is no intersection.Repeat as long as there will be functions for which have no decision on acceptance or discarding.Then the reduction of the amount Bellman functions at each step is determined by the formula: The number of Bellman functions can be reduced according to the expression (17).It was proven, that in the case of signal processing the number of quadratic functions that are required for representation of a Bellman function, generally, does not increase by more than one at each step.Nevertheless, using this approach for processing data based on the two-dimensional the tree approximation of lattice graph neighborhood (Mottl et.al, 1998), the number of quadratic functions on the Bellman functions may be too large and leads to a lack of effective implementation of the procedure.The basic idea of the proposed procedure is to find the groups of closest, in the appropriate sense, quadratic functions using kmeans clustering algorithm, and to replace each of these groups by one quadratic function having the smallest minimun value.We consider quadratic functions in (17) as a points in threedimensional space with coordinates   , , , 1,..., Using k-means clustering algorithm, the distance between quadratic functions can be calculated as the sum of the squares of the difference on the definitive range [a, b]: The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-5/W6, 2015 Photogrammetric techniques for video surveillance, biometrics and biomedicine, 25-27 May 2015, Moscow, Russia Assume that the number k is a predetermined number of groups for k-means clustering algorithm.To preserve the quality of image processing, we do not use directly the final cluster centers which derived by the k-means algorithm.For each of the derived final groups we choose a point that have the lowest third coordinate.
Using the algorithm for reduction in the number of quadratic functions allows you to get the effective implementation of image processing procedure on the basis multi quadratic dynamic programming procedure.Experimental studies show that the vast majority of the original data sets of two or three square functions that are quite fully reflect each function in the Bellman criteria (7).

EXPERIMENTAL RESULTS
In this paper, all the experiments are run on MATLAB 7. MSSIM is defined by the mean of values SSIM: ( ) PSNR (dB) is defined as: SNR (dB) is defined as: and with 3 L  we have following forms of edge functions:

CONCLUSION
Edge preserving MR denoising has become an urgent step in medical imaging to remove noise and to preserve local image features for an accurate and effective diagnosis.In this paper, we proposed a new approach to achieve these aims.The experimental results show that muti-quadratic dynamic programming procedure with the application of the algorithm kmeans for reduction of the number of executed quadratic functions, allows get the high computational efficiency of dynamic programming for image processing, especially in the case of high-resolution images.
14.The test images are 8-bit grayscale brain MR images from Simulated Brain Database (SBD) 1 .In our experiments we compare proposed approaches with the other filters for MR images as Wiener filtering (WF), non -linear total variation (TV), anisotropic diffusion filter (ADF), fast bilateral filter (FBF), the Bayesian least squares with Gaussians Scale Mixture (BLS-GSM).Results are quantified by the Mean Structure Similarity Index (MSSIM), peak-to-signal ratio (PSNR), Signal to Noise Ratio (SNR), and Peak to Signal Noise Ratio (PSNR).
standard deviations (the square root of variance) of images; , P Q  -covariance of images ( , ) P i j 27) On figure 1 show example results of reducing noise MRI with Racian noise level 10   .On figure 2 show example results of reducing noise MRI with Gaussian noise level 15   .

Figure 1 .
Figure 1.Original image (a); Racian noise image (b) with noise level 10   ; denoising results with TV (c), BLS- GSM (d), FBF (e), ADF (f), WF (g) , L=2 (h), L=3 (i) It can be noted that not all of the quadratic function will participate in forming of the final function, because their values are not minimum for any point.Such functions can be dropped using enough simple procedure that takes into account the position of the minimum point and the points of intersection of quadratic functions with each other.