Patent Publication Number: US-8989462-B2

Title: Systems, methods and computer readable storage mediums storing instructions for applying multiscale bilateral filtering to magnetic resonance (RI) images

Description:
CROSS REFERENCE TO RELATED APPLICATIONS 
     This application is a national stage of International Application PCT/US2011/024565, filed Feb. 11, 2011, which claims priority to U.S. Provisional Application No. 61/303,392, filed Feb. 10, 2010, which applications are hereby incorporated by this reference in their entireties. 
     BACKGROUND INFORMATION 
     Magnetic Resonance (MR) images provide invaluable information to medical practitioners. However, due to the presence of imaging artifacts, anatomical variability, varying contrast properties, and poor registration, conventional techniques do not yield satisfactory results over a wide range of scan types and neuroanatomies without manual intervention. Moreover, traditional techniques are not robust enough for large scale analyses. 
     Thus, there is a need for an image processing technique that is robust and accurate. 
     SUMMARY 
     Systems, methods and computer-readable storage mediums storing instructions are provided for segmenting MR images using multiscale bilateral filtering. In some embodiments, a method, performed by a computer having a memory and a processor, for segmenting anatomical landmark in digital medical images, may include: transforming at least one digital image of an anatomical landmark from an Image Domain to a Radon Domain to obtain an original sinogram; decomposing the original sinogram to obtain a plurality of sinograms, each of the plurality of sinograms being at a different scale; combining the plurality of sinograms into a binary sinogram; and reconstructing the binary sinogram into a reconstructed sinogram. 
     In some embodiments, the method may include transforming the reconstructed sinogram from the Radon Domain to the Image Domain to obtain a segmented MR Image. In further embodiments, the method may include displaying the segmented MR Image. In further embodiments, the method may include using the segmented MR Image for attenuation control of a positron emission tomography (PET) Image. In other embodiments, the decomposing is performed using a bilateral filter. 
     In some embodiments, the decomposing may include filtering the original sinogram to a scale of six, wherein the plurality of sinograms is six sinograms. In other embodiments, the method may include filtering each of the plurality of sinograms with a gradient filter after the decomposing the sinogram to the plurality of sinograms. In some embodiments, the gradient filter may be selected based on the anatomical landmark. In other embodiments, the anatomical landmark may include any one, or a combination, of a brain, at least one lung, a thorax, and at least one bone. 
     In further embodiments, the transformation may be based on the following calculation:
 
 P   s (α)= Rf (α, s )=∫ −∞   ∞ ∫ −∞   ∞   f ( x,y )δ( s−x  cos α− y  sin α) dxdy  αε[0 ,π]sεR.  
 
     In some embodiments, a computer-readable storage medium may store instructions for segmenting an anatomical landmark in digital medical images, the instructions may include: transforming at least one digital image of an anatomical landmark from an Image Domain to a Radon Domain to obtain an original sinogram; decomposing the original sinogram to obtain a plurality of sinograms, each of the plurality of sinograms being at a different scale; combining the plurality of sinograms into a binary sinogram; and reconstructing the binary sinogram into a reconstructed sinogram. 
     In some embodiments, the medium may include instructions for transforming the reconstructed sinogram from the Radon Domain to the Image Domain to obtain a segmented MR Image. In further embodiments, the medium may include displaying the segmented MR Image. In further embodiments, the medium may include instructions for using the segmented MR Image for attenuation control of a PET Image. In other embodiments, the decomposing is performed using a bilateral filter. 
     In some embodiments, a system for processing PET and MR images from a PET/MRI machine may include an apparatus including at least one processor; and 
     at least one memory including computer program code. The at least one memory and the computer program code configured to, with the at least one processor, may cause the apparatus to perform at least the following: transforming at least one digital image of an anatomical landmark from an Image Domain to a Radon Domain to obtain an original sinogram; decomposing the original sinogram to obtain a plurality of sinograms, each of plurality of sinograms being at a different scale; combining the plurality of sinograms into a binary sinogram; and reconstructing the binary sinogram into a reconstructed sinogram. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The disclosure can be better understood with the reference to the following drawings and description. The components in the figures are not necessarily to scale, emphasis being placed upon illustrating the principles of the disclosure. 
         FIG. 1  shows a method according to an exemplary embodiment to obtain a segmented tissue image from a scanned original MR image of a subject&#39;s anatomy; 
         FIG. 2  illustrates an example of a schematic diagram of the Radon transform with respect to a brain scan; 
         FIGS. 3(   a ) and  3 ( b ) illustrate an example of a transformation of a MR image of a brain; 
         FIG. 4  illustrates an example of how the Radon transform according to embodiments may improve the signal to noise rate (SNR); 
         FIGS. 5(   a ) and ( b ) illustrate an example of how the Radon transform may decrease the noise in the original images by illustrate the profiles drawing through MR images and corresponding sinograms, respectively; 
         FIG. 6  shows a method of decomposing an original sinogram to a multi-scale sinogram according to an exemplary embodiment; 
         FIG. 7  illustrates an example of a scheme diagram of a bilateral filter; 
         FIGS. 8(   a ) and  8 ( b ) illustrate an example of three scales of generating kernel w i , and its Fourier transforms ŵ i ; 
         FIG. 9  illustrates a scale-space of decomposed images constructed by bilateral filtering; 
         FIG. 10  illustrates a comparison of the profiles of the decomposed images in every scale; 
         FIGS. 11(   a ) and  11 ( b ) show kernels for two filters for a brain scan according to an exemplary embodiment; 
         FIG. 12  illustrates multiscale processing of a brain image to obtain a reconstructed sinogram according to an exemplary embodiment; and 
         FIG. 13  shows an example of an apparatus for providing superposed MR and PET imaging according to an exemplary embodiment. 
     
    
    
     DETAILED DESCRIPTION OF THE EMBODIMENTS 
     The following description, numerous specific details are set forth such as examples of specific components, devices, methods, etc., in order to provide a thorough understanding of embodiments of the disclosure. It will be apparent, however, to one skilled in the art that these specific details need not be employed to practice embodiments of the disclosure. In other instances, well-known materials or methods have not been described in detail in order to avoid unnecessarily obscuring embodiments of the disclosure. While the disclosure is susceptible to various modifications and alternative forms, specific embodiments thereof are shown by way of example in the drawings and will herein be described in detail. It should be understood, however, that there is no intent to limit the disclosure to the particular forms disclosed, but on the contrary, the disclosure is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the disclosure. 
     The methods of the disclosure are described with respect to a brain scan. However, it should be understood that the disclosure is not limited to the brain and may be applied to scans of other anatomical landmarks, including but not limited to, scans of a thorax, at least one lung, and at least one bone. 
     The methods of the disclosure are not limited to the steps described herein. The steps may be individually modified or omitted, as well as additional steps may be added. 
     Unless stated otherwise as apparent from the following discussion, it will be appreciated that terms such as “decomposing,” “filtering,” “combining,” “reconstructing,” “segmenting,” “generating,” “registering,” “determining,” “obtaining,” “processing,” “computing,” “selecting,” “estimating,” “detecting,” “tracking,” “obtaining” or the like may refer to the actions and processes of a computer system, or similar electronic computing device, that manipulates and transforms data represented as physical (e.g., electronic) quantities within the computer system&#39;s registers and memories into other data similarly represented as physical quantities within the computer system memories or registers or other such information storage, transmission or display devices. Embodiments of the methods described herein may be implemented using computer software. If written in a programming language conforming to a recognized standard, sequences of instructions designed to implement the methods may be compiled for execution on a variety of hardware platforms and for interface to a variety of operating systems. In addition, embodiments are not described with reference to any particular programming language. It will be appreciated that a variety of programming languages may be used to implement embodiments of the disclosure. 
     Attenuation correction (AC) is an essential step in the study of quantitative PET imaging of certain body anatomy, specifically those having a well-defined outline such as a scalp, thorax, and lung. One approach is to use MR images for attenuation correction of PET images, involving a segmented tissue technique in which a scan for attenuation correction is taken to segment the scanned anatomy to produce an attenuation correction (AC) map. 
     For example, with respect to a brain scan, the brain may be segmented into the scalp, skull, and brain tissue so that one can assign different attenuation coefficients to the segmented tissue and then produce the AC map. On the AC map, the skull has the largest attenuation coefficient; and thus will dominate the attenuation correction of the image. Once the skull is segmented from MR images, the shape and size of the scalp that is outside of the skull and the brain tissue that is inside of the skull may be determined. 
     Traditionally, many combined PET/MRI machines do not offer transmission scans for attenuation correction because of space limit considerations. Additionally, as mentioned, a number of techniques that have been proposed to segment tissue in MR images fail to provide satisfactory results over a wide range of scan types and anatomies. 
     Segmentation Method 
       FIG. 1  illustrates a method according to an embodiment to obtain a segmented tissue image from a scanned original MR image of a subject&#39;s anatomy. In some embodiments, segmentation method  100  may include a preprocessing step  102  that processes T1-weighted MR images. The preprocessing step  102  may include processing T1-weighted MR images using a predetermined threshold and a Gaussian filter in order to remove the background noise and smooth the original image. The preprocessing may be according to methods well known in the art. 
     Transforming Step 
     The method  100  may further include a transforming step  104  that transforms the MR image from the Image Domain to the Radon Domain to obtain a sinogram. In some embodiments, the Radon transform may be defined as follows from a complete set of line integrals P s (α).
 
 P   s (α)= Rf (α, s )=∫ −∞   ∞ ∫ −∞   ∞   f ( x,y )δ( s−x  cos α− y  sin α) dxdy  αε[0 ,π]sεR.   (1)
 
       FIG. 2  illustrates an example of a schematic diagram of the Radon transform with respect to a brain scan  210 . As shown in  FIG. 2 , s is the perpendicular distance of a line from the origin and α is the angle formed by the distance vector. In this example, the projection image, i.e., sonogram  220  may have two gaps because of low signal for skull. The point ƒ(x, y)  215  in the brain scan  210  corresponds to the sine curve in the Radon domain. According to the Fourier slice theorem, this transformation is invertible. The Fourier slice theorem states that for a 2-D function ƒ(x, y)  215 , the 1-D Fourier transforms of the Radon transform along s, are the 1-D radial samples of the 2-D Fourier transform of ƒ(x, y)  215  at the corresponding angles. This is more fully detailed, for example, in Yang et al., 2010, Proc. SPIE 2010 Mar. 4; 7623, (1):76233K, 7623, 76233K, which is incorporated in its entirety. 
       FIGS. 3(   a ) and  3 ( b ) illustrate an example of a transformation of a MR image  310  of a brain.  FIG. 3(   a ) illustrates the MR image  310  and  FIG. 3(   b ) illustrates the corresponding sinogram  320  in the Radon domain. As shown in  FIG. 3(   b ), there may be two low intensity gaps on top and bottom sides of the sinogram  320  along the vertical direction, which indicate the skull may have relatively low signal intensity on the T1-weighted MR image. 
     Using the Radon Transform according to some embodiments may improve the robustness of the image to noise. Assuming white noise with zero mean is added to the image, the Radon transform of noise will be constant for all of the points and directions and will be equal to the mean value of the noise, which is assumed to be zero, because the Radon transform is line integral of the image, for the continuous case. Accordingly, zero-mean white noise may have no effect on the Radon transform of the image. 
       FIG. 4  illustrates an example of how the Radon transform according to embodiments may improve the signal to noise rate (SNR). Specifically,  FIG. 4  illustrates an example of a brain ellipse model  400  and calculation scheme in the Radon domain. The example may be described by the following equation: 
     
       
         
           
             
               
                 
                   
                     
                       
                         x 
                         2 
                       
                       
                         a 
                         2 
                       
                     
                     + 
                     
                       
                         y 
                         2 
                       
                       
                         b 
                         2 
                       
                     
                   
                   = 
                   
                     1 
                     ⁢ 
                     
                       ( 
                       
                         a 
                         ≥ 
                         b 
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
           
         
       
     
     To demonstrate the improved SNR, the intensity values along the pixels for the Radon transform may be summed as shown in  FIG. 4 . The following calculations are based on the assumption that the point ƒ(x, y) is 2D discrete signals whose intensity values are independent of random variables with mean μ and variance σ 2 . For each point along the projection p s , n s  pixels of ƒ(x, y) are summed and therefore, mean(p s )=n s ·μ and var(p s )=n s ·σ 2 . 
     Meanwhile, n s =s·tan θ=R·sin θ and R=αb/√{square root over (b 2  cos 2  θ+α 2  sin 2  θ)}, where R is the radius of the ellipse area in terms of pixels, and the integer s is the projection index, which varies from −α to α. 
     The average of P s   2  is 
     
       
         
           
             
               
                 
                   
                     A 
                     ⁡ 
                     
                       ( 
                       
                         p 
                         s 
                         2 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       1 
                       
                         2 
                         ⁢ 
                         a 
                       
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         ∑ 
                         
                           s 
                           = 
                           
                             - 
                             a 
                           
                         
                         a 
                       
                       ⁢ 
                       
                         p 
                         s 
                         2 
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     Its expected value is defined as E p =E{A(p s   2 )}. 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           E 
                           p 
                         
                         = 
                         
                           
                             1 
                             
                               2 
                               ⁢ 
                               a 
                             
                           
                           ⁢ 
                           
                             
                               ∑ 
                               
                                 s 
                                 = 
                                 
                                   - 
                                   a 
                                 
                               
                               a 
                             
                             ⁢ 
                             
                               E 
                               ⁢ 
                               
                                 { 
                                 
                                   p 
                                   s 
                                   2 
                                 
                                 } 
                               
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                         
                           
                             1 
                             
                               2 
                               ⁢ 
                               a 
                             
                           
                           ⁢ 
                           
                             
                               ∑ 
                               
                                 s 
                                 = 
                                 
                                   - 
                                   a 
                                 
                               
                               a 
                             
                             ⁢ 
                             
                               [ 
                               
                                 
                                   var 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       p 
                                       s 
                                     
                                     ) 
                                   
                                 
                                 + 
                                 
                                   
                                     ( 
                                     
                                       E 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           p 
                                           s 
                                         
                                         ) 
                                       
                                     
                                     ) 
                                   
                                   2 
                                 
                               
                               ] 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                         
                           
                             1 
                             
                               2 
                               ⁢ 
                               a 
                             
                           
                           ⁢ 
                           
                             
                               ∑ 
                               
                                 s 
                                 = 
                                 
                                   - 
                                   a 
                                 
                               
                               a 
                             
                             ⁢ 
                             
                               ( 
                               
                                 
                                   
                                     n 
                                     s 
                                   
                                   ⁢ 
                                   
                                     σ 
                                     2 
                                   
                                 
                                 + 
                                 
                                   
                                     n 
                                     s 
                                     2 
                                   
                                   ⁢ 
                                   
                                     μ 
                                     2 
                                   
                                 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
     For large α, 
     
       
         
           
             
               
                 
                   
                       
                   
                   ⁢ 
                   
                     
                       
                         E 
                         p 
                       
                       ≈ 
                       
                         
                           1 
                           
                             2 
                             ⁢ 
                             a 
                           
                         
                         ⁢ 
                         
                           
                             ∫ 
                             
                               - 
                               a 
                             
                             a 
                           
                           ⁢ 
                           
                             
                               ( 
                               
                                 
                                   
                                     n 
                                     s 
                                   
                                   ⁢ 
                                   
                                     σ 
                                     2 
                                   
                                 
                                 + 
                                 
                                   
                                     n 
                                     s 
                                     2 
                                   
                                   ⁢ 
                                   
                                     μ 
                                     2 
                                   
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               ⅆ 
                               s 
                             
                           
                         
                       
                     
                     ⁢ 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       ds 
                       = 
                       
                         
                           - 
                           
                             ( 
                             
                               
                                 
                                   a 
                                   3 
                                 
                                 ⁢ 
                                 b 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 sin 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 θ 
                               
                               
                                 
                                   ( 
                                   
                                     
                                       b 
                                       2 
                                     
                                     + 
                                     
                                       
                                         ( 
                                         
                                           
                                             a 
                                             2 
                                           
                                           - 
                                           
                                             b 
                                             2 
                                           
                                         
                                         ) 
                                       
                                       ⁢ 
                                       
                                         sin 
                                         2 
                                       
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       θ 
                                     
                                   
                                   ) 
                                 
                                 
                                   3 
                                   / 
                                   2 
                                 
                               
                             
                             ) 
                           
                         
                         ⁢ 
                         d 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         θ 
                       
                     
                     ⁢ 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     where 
                   
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
             
               
                 
                   
                     E 
                     p 
                   
                   = 
                   
                     
                       ∫ 
                       0 
                       π 
                     
                     ⁢ 
                     
                       
                         
                           ( 
                           
                             
                               
                                 a 
                                 3 
                               
                               ⁢ 
                               b 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               sin 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               θ 
                             
                             
                               
                                 ( 
                                 
                                   
                                     b 
                                     2 
                                   
                                   + 
                                   
                                     
                                       ( 
                                       
                                         
                                           a 
                                           2 
                                         
                                         - 
                                         
                                           b 
                                           2 
                                         
                                       
                                       ) 
                                     
                                     ⁢ 
                                     
                                       sin 
                                       2 
                                     
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     θ 
                                   
                                 
                                 ) 
                               
                               
                                 3 
                                 / 
                                 2 
                               
                             
                           
                           ) 
                         
                         · 
                         
                           ( 
                           
                             
                               
                                 
                                   σ 
                                   2 
                                 
                                 ⁢ 
                                 b 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 sin 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 θ 
                               
                               
                                 
                                   
                                     b 
                                     2 
                                   
                                   + 
                                   
                                     
                                       ( 
                                       
                                         
                                           a 
                                           2 
                                         
                                         - 
                                         
                                           b 
                                           2 
                                         
                                       
                                       ) 
                                     
                                     ⁢ 
                                     
                                       sin 
                                       2 
                                     
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     θ 
                                   
                                 
                               
                             
                             + 
                             
                               
                                 2 
                                 ⁢ 
                                 
                                   μ 
                                   2 
                                 
                                 ⁢ 
                                 
                                   ab 
                                   2 
                                 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   sin 
                                   2 
                                 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 θ 
                               
                               
                                 
                                   
                                     b 
                                     2 
                                   
                                   + 
                                   
                                     
                                       ( 
                                       
                                         
                                           a 
                                           2 
                                         
                                         - 
                                         
                                           b 
                                           2 
                                         
                                       
                                       ) 
                                     
                                     ⁢ 
                                     
                                       sin 
                                       2 
                                     
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     θ 
                                   
                                 
                               
                             
                           
                           ) 
                         
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         ⅆ 
                         θ 
                       
                     
                   
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
           
         
       
     
     In order to simplify this integral, α=b, just like a circular area. 
     
       
         
           
             
               
                 
                   
                     E 
                     p 
                   
                   = 
                   
                     
                       ∫ 
                       0 
                       π 
                     
                     ⁢ 
                     
                       
                         ( 
                         
                           
                             a 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               σ 
                               2 
                             
                             ⁢ 
                             
                               sin 
                               2 
                             
                             ⁢ 
                             θ 
                           
                           + 
                           
                             2 
                             ⁢ 
                             
                               μ 
                               2 
                             
                             ⁢ 
                             
                               θ 
                               2 
                             
                             ⁢ 
                             
                               sin 
                               3 
                             
                             ⁢ 
                             θ 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         ⅆ 
                         θ 
                       
                     
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         
                           E 
                           p 
                         
                         = 
                           
                         ⁢ 
                         
                           
                             a 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               
                                 σ 
                                 2 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     θ 
                                     2 
                                   
                                   - 
                                   
                                     
                                       sin 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       θ 
                                     
                                     4 
                                   
                                 
                                 ) 
                               
                             
                           
                           ⁢ 
                           
                             | 
                             0 
                             π 
                           
                           ⁢ 
                           
                             
                               + 
                               2 
                             
                             ⁢ 
                             
                               μ 
                               2 
                             
                             ⁢ 
                             
                               
                                 a 
                                 2 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     1 
                                     12 
                                   
                                   ⁢ 
                                   
                                     ( 
                                     
                                       
                                         cos 
                                         ⁡ 
                                         
                                           ( 
                                           
                                             3 
                                             ⁢ 
                                             θ 
                                           
                                           ) 
                                         
                                       
                                       - 
                                       
                                         9 
                                         ⁢ 
                                         cos 
                                         ⁢ 
                                         
                                             
                                         
                                         ⁢ 
                                         θ 
                                       
                                     
                                     ) 
                                   
                                 
                                 ) 
                               
                             
                           
                           ⁢ 
                           
                             | 
                             0 
                             π 
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             
                               
                                 π 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 a 
                               
                               2 
                             
                             ⁢ 
                             
                               σ 
                               2 
                             
                           
                           + 
                           
                             
                               8 
                               3 
                             
                             ⁢ 
                             
                               μ 
                               2 
                             
                             ⁢ 
                             
                               a 
                               2 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
           
         
       
     
     Thus for the signal, 
               E   s     ≈           π   ⁢           ⁢   a     2     ⁢     σ   s   2       +       8   3     ⁢     μ   s   2     ⁢       a   2     .               
And for Gaussian white noise with zero mean and variance σ n   2 ,
 
               E   n     ≈         π   ⁢           ⁢   a     2     ⁢       σ   n   2     .             
Thus SNR Radon  in Radon domain may be calculated as
 
     
       
         
           
             
               
                 
                   
                     SNR 
                     Radon 
                   
                   = 
                   
                     
                       
                         E 
                         s 
                       
                       
                         E 
                         n 
                       
                     
                     = 
                     
                       
                         
                           σ 
                           s 
                           2 
                         
                         
                           σ 
                           n 
                           2 
                         
                       
                       + 
                       
                         
                           
                             16 
                             ⁢ 
                             a 
                           
                           
                             3 
                             ⁢ 
                             π 
                           
                         
                         · 
                         
                           
                             μ 
                             s 
                             2 
                           
                           
                             σ 
                             n 
                             2 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
     In the original MR image, the SNR MRI  may be defined as SNR MRI =(σ s   2 +μ s   2 )/σ n   2 . So 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           SNR 
                           Radon 
                         
                         = 
                         
                           
                             
                               
                                 σ 
                                 s 
                                 2 
                               
                               + 
                               
                                 μ 
                                 s 
                                 2 
                               
                             
                             
                               σ 
                               n 
                               2 
                             
                           
                           + 
                           
                             
                               ( 
                               
                                 
                                   
                                     16 
                                     ⁢ 
                                     a 
                                   
                                   
                                     3 
                                     ⁢ 
                                     π 
                                   
                                 
                                 - 
                                 1 
                               
                               ) 
                             
                             · 
                             
                               
                                 μ 
                                 s 
                                 2 
                               
                               
                                 σ 
                                 n 
                                 2 
                               
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                         
                           
                             SNR 
                             MRI 
                           
                           + 
                           
                             
                               ( 
                               
                                 
                                   
                                     16 
                                     ⁢ 
                                     a 
                                   
                                   
                                     3 
                                     ⁢ 
                                     π 
                                   
                                 
                                 - 
                                 1 
                               
                               ) 
                             
                             · 
                             
                               
                                 μ 
                                 s 
                                 2 
                               
                               
                                 σ 
                                 n 
                                 2 
                               
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
     In practice, for real brain MRI, generally, μ s   2 ≧σ s   2 . So 
     
       
         
           
             
               
                 
                   
                     SNR 
                     MRI 
                   
                   = 
                   
                     
                       
                         ( 
                         
                           
                             σ 
                             s 
                             2 
                           
                           + 
                           
                             μ 
                             s 
                             2 
                           
                         
                         ) 
                       
                       / 
                       
                         σ 
                         n 
                         2 
                       
                     
                     ≤ 
                     
                       2 
                       ⁢ 
                       
                         
                           μ 
                           s 
                           2 
                         
                         / 
                         
                           σ 
                           n 
                           2 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       SNR 
                       Radon 
                     
                     ≥ 
                     
                       
                         ( 
                         
                           
                             
                               8 
                               ⁢ 
                               a 
                             
                             
                               3 
                               ⁢ 
                               π 
                             
                           
                           + 
                           
                             1 
                             2 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         SNR 
                         MRI 
                       
                     
                   
                   = 
                   
                     
                       
                         
                           16 
                           ⁢ 
                           a 
                         
                         + 
                         
                           3 
                           ⁢ 
                           π 
                         
                       
                       
                         6 
                         ⁢ 
                         π 
                       
                     
                     ⁢ 
                     
                       SNR 
                       MRI 
                     
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
     Thus, these calculations may demonstrate that SNR Radon  may be at least increased by a factor of 
                   16   ⁢   a     +     3   ⁢   π         6   ⁢   π       .         
If the size of MR image is 256*256, thus
 
               SNR   Radon     ≥         16   ×   256     +     3   ⁢   π         6   ⁢   π       ≈   218.         
So the SNR has been increased by 10 log 10 (218)≈23.4 dB. As a result, this method may be robust to additive noise.
 
       FIGS. 5(   a ) and ( b ) illustrate an example of how the Radon transform may decrease the noise in the original images by illustrate the profiles drawing through MR images and corresponding sinograms, respectively. A first comparison of line profiles  510  is illustrated in  FIG. 5(   a ), in which 8% and 16% noise were added to the original MR image. As shown, the image with 16% noise has severe contamination as the details were lost on the noised image. However, as illustrated in a second comparison of profiles  520  in  FIG. 5(   b ), the profiles between the corresponding sinograms are close, showing that the Radon transform has decreased the noise in the original images. Before the Radon transform, the SNRs are 28.2 dB and 14.4 dB for 8% and 16% noised image, respectively. After the Radon transform, the SNRs are 63.8 dB and 50.1 dB, respectively. Thus, the SNRs are increased by 35.6 dB and 35.7 dB. 
     Decomposing Step 
     The method may further include step  130  that decomposes the sinogram to a multi-scale sinogram. In some embodiments, the decomposing step  130  may performed by the steps shown in  FIG. 6 . The decomposing method may result in a multiscale sinogram representing a series of images with different levels of spatial resolution. In the course scale, the general information is extracted and maintained in images, and in the fine scale, images have more local tissue information. 
     As shown in  FIG. 6 , the decomposing step  130  may include a receiving step  131  for receiving the image in the Radon domain. Because the image is original, the scale, i, equals 0. 
     In some embodiments, the decomposing may be performed by a bilateral filter by a factor of 2. As shown in  FIG. 6 , the decomposing step may further include two filtering steps  132  and  133 . 
     Bilateral filtering is a non-linear filtering technique. Such technique is further detailed in Tomasi and Manduchi, 1998, Computer Vision, 1998, Sixth International Conference, 839-846, which is incorporated by reference in its entirety. In some embodiments, the filter may be a weighted average of the local neighborhood samples, where the weights are computed based on temporal (or spatial in case on images) and radiometric distance between the center sample and the neighboring samples. The bilateral filtering may smooth the images while preserving edges, by means of a nonlinear combination of nearby image values. Bilateral filtering may be described as follows: 
     
       
         
           
             
               
                 
                   
                     
                       h 
                       ⁡ 
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           λ 
                           
                             - 
                             1 
                           
                         
                         ⁡ 
                         
                           ( 
                           x 
                           ) 
                         
                       
                       ⁢ 
                       
                         
                           ∫ 
                           
                             - 
                             ∞ 
                           
                           ∞ 
                         
                         ⁢ 
                         
                           
                             ∫ 
                             
                               - 
                               ∞ 
                             
                             ∞ 
                           
                           ⁢ 
                           
                             
                               I 
                               ⁡ 
                               
                                 ( 
                                 ξ 
                                 ) 
                               
                             
                             ⁢ 
                             
                               
                                 W 
                                 
                                   σ 
                                   s 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   ξ 
                                   - 
                                   x 
                                 
                                 ) 
                               
                             
                             ⁢ 
                             
                               
                                 W 
                                 
                                   σ 
                                   r 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     I 
                                     ⁡ 
                                     
                                       ( 
                                       ξ 
                                       ) 
                                     
                                   
                                   - 
                                   
                                     I 
                                     ⁡ 
                                     
                                       ( 
                                       x 
                                       ) 
                                     
                                   
                                 
                                 ) 
                               
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               ⅆ 
                               ξ 
                             
                           
                         
                       
                     
                   
                   ⁢ 
                   
                       
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
     with the normalization that ensures that the weights for all the pixels add up to one. 
     
       
         
           
             
               
                 
                   
                     
                       λ 
                       ⁡ 
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       
                         ∫ 
                         
                           - 
                           ∞ 
                         
                         ∞ 
                       
                       ⁢ 
                       
                         
                           ∫ 
                           
                             - 
                             ∞ 
                           
                           ∞ 
                         
                         ⁢ 
                         
                           
                             
                               W 
                               
                                 σ 
                                 s 
                               
                             
                             ⁡ 
                             
                               ( 
                               
                                 ξ 
                                 - 
                                 x 
                               
                               ) 
                             
                           
                           ⁢ 
                           
                             
                               W 
                               
                                 σ 
                                 r 
                               
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   I 
                                   ⁡ 
                                   
                                     ( 
                                     ξ 
                                     ) 
                                   
                                 
                                 - 
                                 
                                   I 
                                   ⁡ 
                                   
                                     ( 
                                     x 
                                     ) 
                                   
                                 
                               
                               ) 
                             
                           
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             ⅆ 
                             ξ 
                           
                         
                       
                     
                   
                   ⁢ 
                   
                       
                   
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
           
         
       
     
     Where I(x) and h(x) denote input images and output images, respectively. And W σ     s    measures the geometric closeness between the neighborhood center x and a nearby point ξ and W σ     r    measures the photometric similarity between the pixel at the neighborhood center x and that of a nearby point ξ. Thus, the similarity function W σ     r    operates in the range of the image function I, while the closeness function W σ     s    operates in the domain of I. 
       FIG. 7  illustrates an example of a scheme diagram  700  of a bilateral filter. As shown in the  FIG. 7 , the bilateral filter may replace the pixel value at x with an average of similar and nearby pixel values. In smooth regions, pixel values in a small neighborhood may be similar to each other, and the bilateral filter acts essentially as a standard domain filter, averaging away the small, weakly correlated differences between pixel values caused by noise. This scheme is further described in Elad, 2002, Image Processing, IEEE Transactions, 11(10), 1141-1151, which is incorporated by reference in its entirety. 
     The bilateral filtering may be performed by many different kernels. In some embodiments the bilateral filtering may be performed by a Gaussian filtering. The Gaussian filtering may be shift-invariant Gaussian filtering, in which both the spatial function W σ     s    and the range function W σ     r    are Gaussian functions of the Euclidean distance between their arguments. More specifically, W σ     s    may be described as: 
     
       
         
           
             
               
                 
                   
                     
                       W 
                       
                         σ 
                         s 
                       
                     
                     ⁡ 
                     
                       ( 
                       
                         ξ 
                         - 
                         x 
                       
                       ) 
                     
                   
                   = 
                   
                     ⅇ 
                     
                       - 
                       
                         
                           1 
                           2 
                         
                         ⁡ 
                         
                           [ 
                           
                             
                               d 
                               s 
                               2 
                             
                             / 
                             
                               σ 
                               s 
                               2 
                             
                           
                           ] 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     where d s =∥ξ−x∥ is the Euclidean distance. The range function W σ     r    may be perfectly analogous to W σ     s   : 
     
       
         
           
             
               
                 
                   
                     
                       W 
                       
                         σ 
                         r 
                       
                     
                     ⁡ 
                     
                       ( 
                       
                         ξ 
                         - 
                         x 
                       
                       ) 
                     
                   
                   = 
                   
                     ⅇ 
                     
                       - 
                       
                         
                           1 
                           2 
                         
                         ⁡ 
                         
                           [ 
                           
                             
                               d 
                               r 
                               2 
                             
                             / 
                             
                               σ 
                               r 
                               2 
                             
                           
                           ] 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   16 
                   ) 
                 
               
             
           
         
       
     
     where d r =|I(ξ)−I(x)| is a suitable measure of distance in intensity space. 
     In the scalar case, this may be simply the absolute difference of the pixel difference. The Gaussian range filter may also be insensitive to overall additive changes of image intensity. Thus, the range filter may also be shift-invariant. 
     As discussed above,  FIG. 6  illustrates a multiscale bilateral composition. The decomposing method may further include an obtaining step  135  that obtains each image (sinogram) that is filtered. In further embodiments, the method may further include a determination step  135  that determines whether the filtering step should be repeated so as to build a series of filtered images having different scales as shown in step  136 . The filtering step may be repeated until the number of images at the predetermined scale N is obtained. 
     In some embodiments, a series of filtered images I i  that preserve the strongest edges in I while smoothing small changes in intensity may be obtained based with the following calculations. 
     For the calculations, it is assumed that the original image obtained in step  131  is the 0 th  scale (i=0), that is, set I 0 =I. The image is then filtered through a bilateral filter in steps  132  and  133 . Application of the bilateral filter may compute the following: 
                       I   n     i   +   1       =       1   λ     ⁢       ∑     k   ∈   Ω       ⁢         W       σ   s     ,   i       ⁡     (   k   )       ·       W       σ   r     ,   i       ⁡     (       I     n   +   k     i     -     I   n   i       )       ·     I     n   +   k     i             ⁢     
     ⁢   with   ⁢     
     ⁢     λ   =       ∑     k   ∈   Ω       ⁢         W       σ   s     ,   i       ⁡     (   k   )       ·       W       σ   r     ,   i       ⁡     (       I     n   +   k     i     -     I   n   i       )                     (   17   )               
where n is a pixel coordinate, W σ (x)=exp(−x 2 /σ 2 ), σ s,i  and σ r,i  are the widths of the spatial and range Gaussians respectively and k is an offset relative to n that runs across the support of the spatial Gaussian. The repeated convolution by W σ     s     ,i  may increase the spatial smoothing at each scale i. At the finest scale, the spatial kernel may be set to σ s,i =2 i-1 σ s  (i&gt;0). However, because the bilateral filter is non-linear, the filtered image I i  may not be identical to bilaterally filtering the original input image I with a spatial kernel of cumulative width. The range Gaussian W σ     r     ,i  is an edge-stopping function.
 
     In some embodiments, in order to preserve the edge after several iterations of the bilateral decomposition, σ r,i =2 i-1 σ r  may be set. Increasing the width of the range Gaussian W σ     r     ,i  by a factor of 2 at every scale may increase the chance that an unwanted edge that survives at previous iteration may be smoothed away in later iterations. Also, the initial width σ r  may be set to A/25, where A is the intensity range of the image. 
       FIGS. 8(   a ) and  8 ( b ) illustrate an example of three scales  810 ,  820 , and  830  of generating kernel W i  and its Fourier transforms ŵ i . Dashed lines in  FIG. 8(   a ) show the envelopes of three scales kernel. Dashed lines in  FIG. 8(   b ) indicate smoothing condition boundaries π/2 i-1 &gt;|ω|&gt;π/2 i . 
     As shown in  FIGS. 8(   a ) and  8 ( b ), w i  denotes W σ     s     ,i  and W σ     r     ,i . These generating kernels w i  have a fixed number 2m+1, of non-zero coefficients, yet expand in space by the introduction of 2 i −1 zeros between these coefficients. The expansion of the w i  by inserting zeros leads to the introduction of 2 i  aliases in ŵ i (ω). Thus, the generating kernels themselves do not operate as low-pass filters. However, as long as the following smoothing condition, ŵ i (ω)&lt;ε for π/2 i-1 &gt;|ω|&gt;π/2 i  holds, it will follow that ŵ i (ω)&lt;ε for |ω|&gt;π/2 i  and ŵ i  may be essentially a low-pass filter. 
       FIG. 9  illustrates a scale-space  910  constructed by bilateral filtering. The scale space  910  is composed of a stack of images filtered at different scales where i=0 is the original image.  FIG. 9  illustrates six scales (i=6), however, it would be understood that the decomposing step may be performed to any number of scales. The scale number may be selected based on the anatomy selected to be scanned. 
     As shown in  FIG. 9 , the multiscale bilateral decomposition technique smoothes the images as the scale i increases. The scale i may be considered as the scale level and the original image is at the level 0. When the scale increases, the images become more blurred and contain more general information. Unlike many multi-resolution techniques where the images are down-sampled along the resolution, it is not necessary to subsample the I i  because such downsampling may blur the edges in I i . In addition, downsampling may prevent the decomposition from being translation invariant and could introduce grid artifacts when the coarser scales are manipulated. 
       FIG. 10  illustrates a comparison of the profiles  1000  in every scale. As shown in  FIG. 10 , when the scale increases, small edges in intra-region are smoothed and the big edges in the inter-region are preserved. 
     Gradient Filtering 
     After the sinogram is decomposed into a multiscale sinogram, the segmentation method  100  may further include a filtering step  140  that filters each sinogram through a predetermined gradient filter, as shown in  FIG. 1 . 
     In some embodiments, the images may be filtered through different filters. In some embodiments, the images may be use two sets of filters. In further embodiments, the filters may be selected based on the anatomy scanned.  FIGS. 11(   a ) and  11 ( b ) show kernels  1110  and  1120 , respectively, for two filters utilized for a brain scan. 
     Morphological Image Processing 
     In further embodiments, the method  100  may optionally include a processing step  150  for processing each sinogram with morphological image processing. 
     Combining Filtered Sinograms through Reconstruction of Sinogram 
     In further embodiments, after each sinogram is filtered and processed, the sinograms may be combined in step  160 . The combined sinograms may determine a binary sinogram in step  170 . In some embodiments, a thresholding method is used to determine the binary sinogram from the combined sinograms. The binary sinogram may then be used to reconstruct the sinogram in step  180 . In some embodiments, the reconstruction step  180  may include transforming the sinogram from the Radon Domain to the Image Domain. The transformation may be based on any known calculations. In some embodiments, the reconstruction may be performed using filter back projection or ordered-subset expectation maximization (OSEM) methods. 
       FIG. 12  is an illustration of multiscale processing  1200  performed in steps  160  through  180  on a brain image to obtain a reconstructed sinogram. The scale increases from up to down. After the filtered images are processed, the upper half of the first filtered mage is combined with the lower half of the second filtered image in order to get a new sinogram as shown in ( 1 ) of  FIG. 10 . These images are then used to get a mask in step ( 2 ). From the mask obtained in step ( 2 ), the mask is inversed to obtain the mask in step ( 3 ). After which, the mask from step ( 3 ) may be used to obtain a binary skull in the Radon domain in step ( 4 ). The skull then may be reconstructed in step ( 5 ) to obtain the segmented skull 
     Reconstruction in step  180 , may be performed after obtaining the binary sinogram as discussed above. In some embodiments, a threshold may be utilized in the reconstruction to eliminate some artifacts introduced by reconstruction. Reconstruction may be described as:
 
 f ( x,y )=∫ 0   π   Q   θ ( x  cos θ+ y  sin θ) dθ   (18)
 
     With Q θ  are the ramp-filtered projections. 
     After the reconstruction, in some embodiments, the segmented tissue may be displayed in step  190 . In other embodiments, the segmented MR Image with segmented tissue may be forwarded for further processing. In some embodiments, the segmented MR Image may be used as part of the attenuation control of PET images. In further embodiments, the segmented MR Image with segmented tissue may be both displayed and forwarded. 
     System Implementation 
       FIG. 13  shows an example of a known apparatus  1300  for superposed MR and PET imaging. The apparatus  1300  may include any known MRI tube  102 . The MRI tube  1304  defines a longitudinal direction z, which extends orthogonally to the plane of the drawing of  FIG. 13 . Although not shown, the system may contain a scanner including one or more gamma ray detectors (not shown) (e.g., a ring of gamma ray detectors) incorporated into a RF coil assembly (not shown). The one or more gamma ray detectors may be configured to detect gamma rays from positron annihilations and may include a plurality of scintillators and photodetectors arranged circumferentially about a gantry (i.e., a ring of gamma ray detectors or a detector ring). 
     As shown in  FIG. 13 , the system may further include a plurality of PET detection units  1306  arranged in mutually opposing pairs around the longitudinal direction z are arranged coaxially within the MRI tube  1304 . The PET detection units  1306  may preferably include a photodiode array  1310 , such an APD photodiode array, with an upstream array of crystals  1308 , such as LSO crystals, and an electrical amplifier circuit (AMP)  1306 . 
     The apparatus may further include a computer system  1320  to carry out the image processing for superposed MR and PET imaging. The computer system  1320  may further be used to control the operation of the system or a separate system may be included. 
     The computer system  1320  may include a number of modules that communicate with each other through electrical and/or data connections (not shown). Data connections may be direct wired links or may be fiber optic connections or wireless communications links or the like. The computer system  1320  may also be connected to permanent or back-up memory storage, a network, or may communicate with a separate system control through a link (not shown). The modules may include a CPU  1322 , a memory  1324 , an image processor  1326 , an input device  1328 , a display  1330 , and a printer interface  1332 . The computer system  1320  may also be connected to another computer system as well as a network. 
     The CPU  1322  may be one or more of any known central processing unit, including but not limited to a processor, or a microprocessor. The CPU  1322  may be coupled directly or indirectly to memory elements. The memory  1324  may include random access memory (RAM), read only memory (ROM), disk drive, tape drive, etc., or a combinations thereof). The memory may also include a frame buffer for storing image data arrays. 
     The described processes (e.g.,  FIGS. 1 and 6 ) may be implemented as a routine that is stored in memory  1324  and executed by the CPU  1322 . As such, the computer system  1320  may be a general purpose computer system that becomes a specific purpose computer system when executing the routine of the disclosure. The computer system  1320  may also include an operating system and micro instruction code. The various processes and functions described herein may either be part of the micro instruction code or part of the application program or routine (or combination thereof) that is executed via the operating system. In addition, various other peripheral devices may be connected to the computer platform such as an additional data storage device, a printing device. and I/O devices. 
     The image processor  1326  may be one or more of any known central processing unit, including but not limited to a processor, or a microprocessor. In some embodiments, the image processor  1326  also processes the data. In other embodiments, the image processor  1326  may be replaced by image processing functionality on the CPU  1322 . 
     The input device  1328  may include a mouse, joystick, keyboard, track ball, touch activated screen, light wand, voice control, or any similar or equivalent input device, and may be used for interactive geometry prescription. The input device  1328  may control the production, display of images on the display  1330 , and printing of the images by the printer interface  1332 . The display  1330  and the printer interface  1332  may be any known display screen and the printer interface  1332  may be any known printer, either locally or network connected. 
     An MRI scan is complete when one or more sets of raw k-space data has been acquired in by the image processor  1326 . The image processor  1326  reconstructs the raw k-space data by transforming the data raw k-space data (via Fourier transformation or another technique) into image data. This image data may then be stored in the memory  1324 . In other embodiments, another computer system may assume the duties of image reconstruction or other functions of the image processor  1326 . In response to commands received from the input device  1328 , the image data stored in the memory  1324  may be archived in long term storage or may be further processed by the image processor  1326  and presented on the display  1330 . PET images may be reconstructed by the image processor  1326  and may be combined with MR images to produce hybrid structural and metabolic or functional images. 
     It is to be understood that the embodiments of the disclosure be implemented in various forms of hardware, software, firmware, special purpose processes, or a combination thereof. In one embodiment, the disclosure may be implemented in software as an application program tangible embodied on a computer readable program storage device. The application program may be uploaded to, and executed by, a machine comprising any suitable architecture. The system and method of the present disclosure may be implemented in the form of a software application running on a computer system, for example, a mainframe, personal computer (PC), handheld computer, server, etc. The software application may be stored on a recording media locally accessible by the computer system and accessible via a hard wired or wireless connection to a network, for example, a local area network, or the Internet. 
     It is to be further understood that, because some of the constituent system components and method steps depicted in the accompanying figures can be implemented in software, the actual connections between the systems components (or the process steps) may differ depending upon the manner in which the disclosure is programmed. Given the teachings of the disclosure provided herein, one of ordinary skill in the related art will be able to contemplate these and similar implementations or configurations of the disclosure. 
     While the disclosure has been described in detail with reference to exemplary embodiments, those skilled in the art will appreciate that various modifications and substitutions can be made thereto without departing from the spirit and scope of the disclosure as set forth in the appended claims. For example, elements and/or features of different exemplary embodiments may be combined with each other and/or substituted for each other within the scope of this disclosure and appended claims.