Abstract:
Systems and methods for sonar image feature extraction, which fit a parameterized statistical model to sonar image data to extract features that describe image textures for follow-on pattern classification tasks. The systems and methods estimate the parameters of an assumed statistical model from the pixel values in the sonar images. The parameters of the distribution can then be used as features in a discriminant classifier to perform tasks such as texture segmentation or environmental characterization. The systems and methods divide a sonar image into separate overlapping cells. The ACF parameters of each cell are sequentially estimated using various methods. Based on the ACF parameters obtained, each cell is assigned a label via a pattern classification algorithm. The labeled image can be used in target recognition, environmental characterization and texture segmentation.

Description:
STATEMENT OF GOVERNMENT INTEREST 
     The invention described herein may be manufactured and used by or for the Government of the United States of America for governmental purposes without the payment of any royalties thereon or therefore. 
    
    
     CROSS REFERENCE TO OTHER PATENT APPLICATIONS 
     Not applicable. 
     BACKGROUND OF THE INVENTION 
     (1) Field of the Invention 
     The present invention relates to sonar image feature extraction and more specifically to fitting a parameterized statistical model to sonar image data to extract features that describe image textures for follow-on pattern classification tasks. 
     (2) Description of the Prior Art 
     As is known in the art, parameterized statistical modeling of sonar imagery can produce simple, robust quantitative descriptors for complex textured seabed environments. One such known model is a correlated K-distribution model for synthetic aperture radar (SAR) image pixels, with parameters defining spatial correlation lengths in vertical and horizontal dimensions along with a spatial frequency component in one direction to model periodic variation in the autocorrelation function (ACF). 
     The K-distribution has been shown to adequately describe the single-point amplitude and intensity pixel statistics of sonar images. One known method to derive the K-distribution probability density function (pdf) equation is to treat the sonar intensity scene as being composed of exponentially-distributed speckle modulated by a gamma-distributed parameter in the form of a product of two random variables. Along with its intuitive appeal, this convention allows for the introduction of correlated texture through a correlated gamma distribution. 
     As is known to those of skill in the art, the analytic expression for the pdf for a K-distributed sonar image intensity pixel I(x,y) is: 
                           p   I     ⁡     (   x   )       =       α     Γ   ⁡     (   v   )         ⁢       (     α   2     )     v     ⁢       (   x   )         v   -   1     2       ⁢       K     v   -   1       ⁡     (     α   ⁢     x       )       ⁢     u   ⁡     (   x   )           ;     v   &gt;   0       ,           (   1   )               
where α is the scale parameter, ν is the shape parameter of the distribution and K ν-1  is the modified Bessel function of the second kind.
 
     The pixel statistical parameter estimation described above relies on independence between pixel samples. Thus, the correlation between neighboring pixels is ignored. As such, the pdf in Equation (1) describes the one-dimensional statistics of the image pixel values, but does not describe texture. 
     What is needed is a parameterized statistical model that can be fitted to synthetic aperture sonar image data to extract features or descriptors for follow-on pattern classification tasks. What is further needed is for the model to extract parameters that describe image textures such as correlation length, periodicity or ripple, orientation, and other statistical descriptors for use in automatic target recognition, through-the-sensor environmental characterization, and texture segmentation schemes. 
     SUMMARY OF THE INVENTION 
     It is therefore a general purpose and primary object of the present invention to provide a sonar image feature extractor wherein a parameterized statistical model is fitted to sonar image data to extract image texture features, such as correlation length, periodicity or ripple, orientation and other statistical descriptors. Such features can be used in automatic target recognition, through-the-sensor environmental characterization and texture segmentation schemes. 
     The sonar image texture segmentation model described herein estimates the parameters of an assumed statistical model from the pixel values in the sonar images. The parameters of the distribution can then be used as features in a discriminant classifier to perform tasks such as texture segmentation or environmental characterization. 
     The model expands the known correlated K-distribution model into a Gaussian mixture model to allow more degrees of freedom including modeling periodicity. The model adds an angular component to parameterize rotation of the autocorrelation function (ACF) about the vertical and horizontal axes. The model further provides a method to estimate texture parameters using an Expectation-Maximization (EM) method for truncated data. 
     As will be explained in further detail hereinafter, implementing the model includes dividing a sonar image into separate overlapping cells. The overlap amount is a calculated tradeoff between texture parameter resolution and processing speed. The ACF parameters are sequentially estimated using various methods. Based on the ACF parameters obtained, each cell is assigned a label via a pattern classification algorithm. The labeled image can be used in target recognition, environmental characterization and texture segmentation. 
     The parameters of the ACF of the imaging point spread function (PSF) are estimated using least squares. The shape parameter of the single-point image pixel pdf is estimated using a method known in the art for estimating parameters of K-distributed clutter. The remaining parameters are estimated using the EM method. 
     In one embodiment, a method of texture segmentation for a sonar image, S(x,y), includes obtaining initial estimates of a first set of parameters of the sonar image, obtaining an estimate of an intensity autocorrelation function (ACF) of the image, and initializing a second set of parameters, Θ i ={α i , [μ X     i    μ Y     i   ], Σ i }, of the sonar image based on the first set of parameters. Additionally, the method includes initializing a likelihood function based on the first and second sets of parameters, maximizing the likelihood function to obtain updated parameters, calculating an updated intensity ACF based on the updated parameters, and iteratively returning to maximizing the likelihood function until convergence of the updated intensity ACF. 
     The step of obtaining the estimate of the intensity ACF includes estimating the intensity ACF via discrete two-dimensional Fourier transform: 
                   R   ^     I     ⁡     (     X   ,   Y     )       =         1   NO     ·     F     -   1         ⁢     {     F   ⁢     {                S   ⁡     (     x   ,   y     )            2     ⁢     *   •     ⁢     conj   (     F   ⁢     {            S   ⁡     (     x   ,   y     )       ⁢     ❘   2       }     )       }       ,                   
for a (N×O) cell of the image. Also, the step of obtaining initial estimates includes estimating a mean intensity,
 
                   μ   ^     I     =       1   NO     ⁢       ∑     x   i       ⁢       ∑     y   i       ⁢     S   ⁡     (       x   i     ,     y   i       )               ,         
for the (N×O) cell, and estimating an imaging point spread function (PSF) ACF, {circumflex over (R)} h (X, Y), for the cell via discrete two-dimensional Fourier transform based on said mean intensity:
 
                   R   ^     h     ⁡     (     X   ,   Y     )       =             (       μ   ^     I     )       -   1       NO     ·     F     -   1         ⁢       {     F   ⁢     {     S   ⁡     (     x   ,   y     )       }     ⁢     *   •     ⁢     conj   ⁡     (     F   ⁢     {     S   ⁡     (     x   ,   y     )       }       )         }     .             
Obtaining the initial estimates further includes estimating PSF definition variables, {circumflex over (β)} x  and {circumflex over (β)} y , via least squares matrix computation based on {circumflex over (R)} h (X, Y) and {circumflex over (μ)} I , and estimating shape parameter, {circumflex over (ν)}, via K-distributed clutter parameter estimation.
 
     The step of maximizing the likelihood function includes calculating a first likelihood parameter E j {τ i   (j) ([  x j   ,  y j   ]|Θ i   (j-1) )} from the expression 
                   〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j     =         α   i     (     k   -   1     )       ⁢       G   i     ⁡     (         [         x   j     _     ,       y   j     _       ]     |     [       μ     X   i       (     k   -   1     )       ⁢     μ     Y   i       (     k   -   1     )         ]       ,     Σ   i     (     k   -   1     )         )             Σ     i   =   1     M     ⁢     α   i     (     k   -   1     )       ⁢       G   i     ⁡     (         [         x   j     _     ,       y   j     _       ]     |     [       μ     X   i       (     k   -   1     )       ⁢     μ     Y   i       (     k   -   1     )         ]       ,     Σ   i     (     k   -   1     )         )             ,         
where
 
                 G   i     ⁡     (       [       μ     X   i       ⁢     μ     Y   i         ]     ,     Σ   i       )       =         (     2   ⁢           ⁢   π     )       -   1       ⁢            Σ   i                -       1   2     ⁡     [       [   XY   ]     -     [       μ     x   i       ⁢     μ     y   i         ]       ]         ⁢       Σ   i     -   1       ⁡     [         [   XY   ]     T     -       [       μ     x   i       ⁢     μ     y   i         ]     T       ]         ,               
calculating a scale parameter
 
                 α   i     (   k   )       =         Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j           Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )             ,         
where
 
                 m   j     ⁡     (     Θ     (   k   )       )       =     {             n   j     ,     j   =   1     ,   …   ⁢           ,   r                 n   ·         P   j     ⁡     (     Θ   j     (   k   )       )         P   ⁡     (     Θ     (   k   )       )           ,     j   =     r   +   1       ,   …   ⁢           ,   v   ,     and   ⁢           ⁢         P   j     ⁡     (   Θ   )         P   ⁡     (   Θ   )                           
is a bin probability of an area defined by a resolution cell (x j ,y j ), calculating a vector of component means
 
                   [       μ     X   i       (   k   )       ⁢     μ     Y   i       (   k   )         ]     T     =         Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢           〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j     ⁡     [         x   j     _     ,       y   j     _       ]       T           Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j           ,         
and calculating a matrix variable
 
                 Σ   i     (   k   )       =         Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j     ⁢     K   T     ⁢   K         Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j           ,         
where K=[[  x j   ,  y j   ]−[μ X     i     (k)  μ Y     i     (k) ]] and Σ i  comprises correlation length parameters l x     i    and l y     i    and ACF rotation parameter θ i . Maximizing the likelihood function further includes obtaining an updated likelihood function based on the first likelihood parameter, the scale parameter, the vector of component means, and the matrix variable. Further, the method iteratively returns to calculating the first likelihood parameter, the scale parameter, the vector of component means and the matrix variable until convergence of said updated likelihood function. The updated parameters are obtained based on the updated likelihood function.
 
     Additionally, the cell can be a subset of the image and the method can further include dividing the image into a plurality of cells, wherein the cells overlap. Cell parameters associated with the updated intensity ACF for each of the plurality of cells can be stored and textures of the image based on the stored cell parameters can be classified. 
     In another embodiment, a method of segmenting textures in a sonar image, S(x,y), includes dividing the image into a plurality of (N×O) cells, wherein the cells overlap, and obtaining initial estimates of a first set of parameters of a first one of the cells. The method also includes obtaining an estimate of an intensity autocorrelation function (ACF), {circumflex over (R)} I (X, Y), of the first cell via discrete two-dimensional Fourier transform: 
     
       
         
           
             
               
                 
                   R 
                   ^ 
                 
                 I 
               
               ⁡ 
               
                 ( 
                 
                   X 
                   , 
                   Y 
                 
                 ) 
               
             
             = 
             
               
                 
                   1 
                   NO 
                 
                 · 
                 
                   F 
                   
                     - 
                     1 
                   
                 
               
               ⁢ 
               
                 { 
                 
                   F 
                   ⁢ 
                   
                     { 
                     
                       
                         
                            
                           
                             S 
                             ⁡ 
                             
                               ( 
                               
                                 x 
                                 , 
                                 y 
                               
                               ) 
                             
                           
                            
                         
                         2 
                       
                       ⁢ 
                       
                         * 
                         • 
                       
                       ⁢ 
                       
                         
                           conj 
                           ( 
                           
                             F 
                             ⁢ 
                             
                               { 
                               
                                  
                                 
                                   
                                     S 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         x 
                                         , 
                                         y 
                                       
                                       ) 
                                     
                                   
                                   ⁢ 
                                   
                                     ❘ 
                                     2 
                                   
                                 
                                 } 
                               
                               ) 
                             
                           
                           } 
                         
                         . 
                       
                     
                   
                 
               
             
           
         
       
     
     Additionally, the method includes initializing a second set of parameters, Θ i ={α i , [μ X     i    μ Y     i   ], Σ i }, of the first cell based on the first set of parameters, initializing a likelihood function based on the first and second sets of parameters, and maximizing the likelihood function to obtain updated parameters. Further, the method calculates an updated intensity ACF based on the updated parameters, iteratively returns to maximizing the likelihood function until convergence of the updated intensity ACF, and stores the estimated ACF parameters in a matrix. The method iteratively returns to obtaining initial estimates of a first set of parameters for each of the plurality of cells and segments textures of the sonar image based on the estimated texture parameters of the cells. 
     The step of obtaining initial estimates further includes estimating a mean intensity, 
                 μ   I     =       1   NO     ⁢       ∑     x   i       ⁢       ∑     y   i       ⁢     S   ⁡     (       x   i     ,     y   i       )               ,         
and estimating an imaging point spread function (PSF) ACF, {circumflex over (R)} h (X, Y), for the cell via discrete two-dimensional Fourier transform based on said mean intensity:
 
                   R   ^     h     ⁡     (     X   ,   Y     )       =             (       μ   ^     I     )       -   1       NO     ·     F     -   1         ⁢       {     F   ⁢     {     S   ⁡     (     x   ,   y     )       }     ⁢     *   •     ⁢     conj   ⁡     (     F   ⁢     {     S   ⁡     (     x   ,   y     )       }       )         }     .             
The step of obtaining initial estimates also includes estimating PSF definition variables, {circumflex over (β)} x  and {circumflex over (β)} y , via least squares matrix computation based on {circumflex over (R)} h (X, Y) and {circumflex over (μ)} I , and estimating shape parameter, {circumflex over (ν)}, via K-distributed clutter parameter estimation.
 
     The step of maximizing the likelihood function further includes calculating likelihood parameters, obtaining an updated likelihood function based on the likelihood parameters, iteratively returning to calculating the likelihood parameters until convergence of the updated likelihood function, and obtaining the updated parameters based on the updated likelihood function. The step of calculating likelihood parameters further includes calculating a first likelihood parameter E j {τ i   (j) ([  x j   ,  y j   ]|Θ i   (j-1) )} from the expression 
                   〈       τ   i     (   k   )       ⁢     |     ⁢     Θ   i     (     k   -   1     )         〉     j     =         α   i     (     k   -   1     )       ⁢       G   i     ⁡     (         [         x   j     _     ,       y   j     _       ]     ⁢       |     ⁡     [       μ     X   i       (     k   -   1     )       ⁢     μ     Y   i       (     k   -   1     )         ]         ,     ∑   i     (     k   -   1     )         )             ∑     i   =   1     M     ⁢       α   i     (     k   -   1     )       ⁢       G   i     ⁡     (         [         x   j     _     ,       y   j     _       ]     ⁢       |     ⁡     [       μ     X   i       (     k   -   1     )       ⁢     μ     Y   i       (     k   -   1     )         ]         ,     ∑   i     (     k   -   1     )         )               ,         
where
 
                 G   i     (       [       μ     X   i       ⁢     μ     Y   i         ]     ,     ∑   i       )     =         (     2   ⁢   π     )       -   1       ⁢            ∑   i            -     1   2         ⁢       ⅇ       -       1   2     ⁡     [       [     X   ⁢           ⁢   Y     ]     -     [       μ     x   i       ⁢     μ     y   i         ]       ]         ⁢       ∑   i     -   1       ⁢     [         [     X   ⁢           ⁢   Y     ]     T     -       [     ⁣       μ     x   i       ⁢     μ     y   i           ]     T       ]           .             
Also, calculating the likelihood parameters includes calculating a scale parameter
 
                 α   i     (   k   )       =         ∑     j   =   1     v     ⁢         m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       ⁢     |     ⁢     Θ   i     (     k   -   1     )         〉     j             ∑     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )             ,         
where
 
                 m   j     ⁡     (     Θ     (   k   )       )       =     {             n   j     ,     j   =   1     ,   …   ⁢           ,   r                 n   ·         P   j     ⁡     (     Θ   j     (   k   )       )         P   ⁡     (     Θ     (   k   )       )           ,     j   =     r   +   1       ,   …   ⁢           ,   v   ,     and   ⁢           ⁢         P   j     ⁡     (   Θ   )         P   ⁡     (   Θ   )                           
is a bin probability of an area defined by a resolution cell (x j ,y j ).
 
     The step of calculating the likelihood parameters further includes calculating a vector of component means 
                   [       μ     X   i       (   k   )       ⁢     μ     Y   i       (   k   )         ]     T     =         ∑     j   =   1     v     ⁢         m   j     ⁡     (     Θ     (   k   )       )       ⁢           〈       τ   i     (   k   )       ⁢     |     ⁢     Θ   i     (     k   -   1     )         〉     j     ⁡     [         x   j     _     ,       y   j     _       ]       T             ∑     j   =   1     v     ⁢         m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       ⁢     |     ⁢     Θ   i     (     k   -   1     )         〉     j             ,         
and calculating a matrix variable
 
                 ∑   i     (   k   )       ⁢     =         ∑     j   =   1     v     ⁢         m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       ⁢     |     ⁢     Θ   i     (     k   -   1     )         〉     j     ⁢     K   T     ⁢   K           ∑     j   =   1     v     ⁢         m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       ⁢     |     ⁢     Θ   i     (     k   -   1     )         〉     j               ,         
where K=[[  x j   ,  y j   ]−[μ X     i     (k)  μ Y     i     (k) ]] and Σ i  comprises correlation length parameters, l x     i    and l y     i   , and ACF rotation parameter θ i .
 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       A more complete understanding of the invention and many of the attendant advantages thereto will be readily appreciated as the same becomes better understood by reference to the following detailed description when considered in conjunction with the accompanying drawings wherein like references numerals and symbols designate identical or corresponding parts throughout the several views and wherein: 
         FIG. 1  illustrates a block diagram of a method for implementing sonar image feature segmentation; 
         FIG. 2  illustrates a block diagram of a parameter estimation portion of the method of  FIG. 1 ; 
         FIG. 3  illustrates a block diagram of an initial estimation portion of the method of  FIG. 1 ; 
         FIG. 4  illustrates a block diagram a likelihood function calculating portion of the method of  FIG. 1 ; 
         FIGS. 5A-D  illustrate synthetic aperture sonar image textures; 
         FIGS. 6A-D  illustrate log intensity ACF estimates associated with the textures depicted in  FIGS. 5A-5D , respectively; 
         FIGS. 7A-D  illustrate reconstructed log intensity ACFs resulting from applying the method of  FIG. 1  to the textures depicted in  FIGS. 5A-5D , respectively; 
         FIG. 8A  illustrates an exemplary image having two synthetic textures; and 
         FIGS. 8B and 8C  illustrate the respective {circumflex over (l)} x  and {circumflex over (l)} y  parameter estimations for the image of  FIG. 8A . 
     
    
    
     DESCRIPTION OF THE INVENTION 
     As previously described herein with relation to Equation (1), the compound K-distribution has been used to model sonar image pixel statistics. In deriving Equation (1), the sonar intensity scene is treated as being composed of exponentially-distributed speckle modulated by a gamma-distributed parameter in the form of a product of two random variables. Along with its intuitive appeal, this convention allows for the introduction of correlated texture through a correlated gamma distribution. However, as also previously noted, pixel statistical parameter estimation in Equation (1) assumes independence between pixel samples and correlation between neighboring pixels is ignored. 
     To describe two-dimensional correlations in the imagery, the following known expression for the normalized intensity ACF can be used: 
                       R     I   N       ⁡     (     X   ,   Y     )       =             R   I     ⁡     (     X   ,   Y     )         μ   I   2       -   1     =         [       R   h     ⁡     (     X   ,   Y     )       ]     2     +       1   v     ·             ∑   i     ⁢         η   i     ·     l     x   i         ⁢           ⁢     l     y   i       ⁢           ⁢              ∑   i            1   2       ·     [           ⁢               ⅇ         1   2     ⁡     [       [     X   ⁢           ⁢   Y     ]     -     [       μ     X   i       ⁢     μ     Y   i         ]       ]       ⁢       ∑   i     -   1       ⁢     [         [     X   ⁢           ⁢   Y     ]     T     -       [       μ     X   i       ⁢     μ     Y   i         ]     T       ]           +       [       R   h     ⁡     (     X   ,   Y     )       ]     2       ]     ,                             (   2   )               
where Σ i =Ψ i Λ i Ψ i   T . In this case, Λ i  is a diagonal matrix of mixture component correlation length parameters
 
                       Λ   i     =     [           l     x   i     2         0           0         l     y   i     2           ]       ,           (   3   )               
and ψ i  is a mixture component rotation matrix
 
                       Ψ   i     =     [           cos   ⁢           ⁢     θ   i             sin   ⁢           ⁢     θ   i                   -   sin     ⁢           ⁢     θ   i             cos   ⁢           ⁢     θ   i             ]       ,           (   4   )               
where θ i  is the counterclockwise rotation of the sonar cross section (SCS) ACF mixture component from the X-axis. To summarize, Σ i  comprises correlation length parameters l x     i    and l y     i    and ACF rotation parameter θ i .
 
     Other parameters of interest in Equation (2) are as follows: 
     μ X     i    and μ Y     i    are vectors of component means in the X-spatial and Y-spatial lag directions, respectively. In rippled textures, μ X     i    and μ Y     i    values correspond to the spatial period of the rippled texture; 
     R h (X, Y) is the autocorrelation function of the imaging point spread function. The variables fix β x  and β y  define this function, where 
                 R   h     ⁡     (     X   ,   Y     )       ⁢     =   Δ     ⁢         ⅇ         -     1   8       ·     [     X   ⁢           ⁢   Y     ]       ⁢         B     -   1       ⁡     [     X   ⁢           ⁢   Y     ]       T         ⁢           ⁢   and   ⁢           ⁢   B     =       [           β   x   2         0           0         β   y   2           ]     .             
The parameters β x  and β y  are directly proportional to the synthetic aperture sonar (SAS) image pixel resolution;
 
     μ I  is the mean of the intensity image. The intensity ACF is normalized by this quantity so textures in images of different scales can be compared equally; 
     η i  is the mixture component proportion; and 
     ν is a shape parameter of single-point image pixel pdf. 
     As will be described in further detail hereinafter, the parameters of Equation (2) are estimated by sequentially estimating the parameters of R h (X, Y) using least-squares, ν using a known method for estimating parameters of K-distributed clutter, and the remaining parameters using an Expectation Maximization algorithm. Given estimates of the intensity ACF {circumflex over (R)} I , the mean intensity {circumflex over (μ)} I , and the estimated parameters {circumflex over (β)} x , {circumflex over (β)} y , and {circumflex over (ν)}, the remaining parameters l x     i   , l y     i   , [μ X     i   , μ Y     i   ], θ i , η i , i={1, 2, . . . , M} can be found by manipulating Equation (2), where M is the number of Gaussian mixture components. 
     Dropping reference to spatial lag (X, Y) to provide for compact notation hereafter, manipulation of Equation (2) yields 
                           R   ^       I   N       (     k   +   1     )       -         R   ^     h   2     ⁡     [     1   +       i   v     ⁢       f     (   k   )       ⁡     (       η   i     ,     l     x   i       ,     l     y   i         )           ]         =       1     v   ^       ⁢         f     (     k   +   1     )       ⁡     (       η   i     ,     l     x   i       ,     l     y   i         )       ·     [     ⅇ         1   2     ⁡     [       [     X   ⁢           ⁢   Y     ]     -     [       μ     X   i       ⁢     μ     Y   i         ]       ]       ⁢         (     ∑     (     k   +   1     )       )       -   1       ⁡     [         [     X   ⁢           ⁢   Y     ]     T     -       [       μ     X   i       ⁢     μ     Y   i         ]     T       ]           ]           ,           (   5   )               
where
 
               f   ⁡     (       η   i     ,     l     x   i       ,     l     y   i         )       ⁢     =   Δ     ⁢       ∑   i     ⁢         η   i     ·     l     x   i         ⁢     l     y   i       ⁢              ∑   i            1   2       .               
Those of skill in the art will recognize that the right hand side of Equation (5) is a scaled M-component Gaussian mixture model of the form
 
                     C   ·       ∑     i   =   1     M     ⁢       α   i     ⁢       G   i     (       [       μ     X   i       ⁢     μ     Y   i         ]     ,     ∑   i       )           ,   where           (   6   )               C   =     2   ⁢     π   ·         ∑   j     ⁢       η   j     ⁢     l     x   j       ⁢     l     y   j             v   ^         ⁢           ⁢   and             (   7   )                   α   i     =         η   i     ·     l     x   i       ·     l     y   i             ∑   j     ⁢       η   j     ⁢     l     x   j       ⁢     l     y   j               ,           (   8   )               
so that Σ i α i =1 and
 
     
       
         
           
             
               
                 G 
                 i 
               
               ⁡ 
               
                 ( 
                 
                   
                     [ 
                     
                       
                         μ 
                         
                           X 
                           i 
                         
                       
                       ⁢ 
                       
                         μ 
                         
                           Y 
                           i 
                         
                       
                     
                     ] 
                   
                   , 
                   
                     ∑ 
                     i 
                   
                 
                 ) 
               
             
             = 
             
               
                 
                   ( 
                   
                     2 
                     ⁢ 
                     π 
                   
                   ) 
                 
                 
                   - 
                   1 
                 
               
               ⁢ 
               
                 
                    
                   
                     ∑ 
                     i 
                   
                    
                 
                 
                   1 
                   2 
                 
               
               ⁢ 
               
                 
                   ⅇ 
                   
                     
                       
                         1 
                         2 
                       
                       ⁡ 
                       
                         [ 
                         
                           
                             [ 
                             
                               X 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               Y 
                             
                             ] 
                           
                           - 
                           
                             [ 
                             
                               
                                 μ 
                                 
                                   x 
                                   i 
                                 
                               
                               ⁢ 
                               
                                 μ 
                                 
                                   y 
                                   i 
                                 
                               
                             
                             ] 
                           
                         
                         ] 
                       
                     
                     ⁢ 
                     
                       
                         ∑ 
                         i 
                         
                           - 
                           1 
                         
                       
                       ⁢ 
                       
                         [ 
                         
                           
                             
                               [ 
                               
                                 X 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 Y 
                               
                               ] 
                             
                             T 
                           
                           - 
                           
                             
                               [ 
                               
                                 
                                   μ 
                                   
                                     x 
                                     i 
                                   
                                 
                                 ⁢ 
                                 
                                   μ 
                                   
                                     y 
                                     i 
                                   
                                 
                               
                               ] 
                             
                             T 
                           
                         
                         ] 
                       
                     
                   
                 
                 . 
               
             
           
         
       
     
     Using the known EM algorithm for truncated two-dimensional Gaussian mixture models, the parameters Θ (k) ={α 1   (k) , . . . , α M   (k) , └μ X     1     (k)  μ Y     1     (k) ┘, . . . , └μ X     M     (k)  μ Y     1     (k) ┘, Θ 1   (k) , . . . , Σ M   (k) } at step k are found by maximizing the likelihood function 
                     L   ⁡     (     Θ     (   k   )       )       =       ∑   X     ⁢       ∑   Y     ⁢           g     (   k   )       ⁡     (     X   ,   Y     )       ·   ln     ⁢           ⁢       C     (   k   )       ·       ∑     i   =   1     M     ⁢       α   i     (   k   )       ⁢       G   i     (   k   )       ⁡     (       [       μ     X   i       (   k   )       ⁢     μ     Y   i       (   k   )         ]     ,     ∑   i     (   k   )         )                           (   9   )               
where
 
     
       
         
           
             
               
                 
                   
                     
                       g 
                       
                         ( 
                         k 
                         ) 
                       
                     
                     ⁡ 
                     
                       ( 
                       
                         X 
                         , 
                         Y 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               R 
                               ^ 
                             
                             I 
                           
                           
                             
                               μ 
                               ^ 
                             
                             σ 
                             2 
                           
                         
                         - 
                         1 
                         - 
                         
                           
                             
                               R 
                               ^ 
                             
                             h 
                             2 
                           
                           ⁡ 
                           
                             [ 
                             
                               1 
                               + 
                               
                                 
                                   1 
                                   
                                     v 
                                     ^ 
                                   
                                 
                                 ⁢ 
                                 
                                   f 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       
                                         η 
                                         i 
                                         
                                           ( 
                                           
                                             k 
                                             - 
                                             1 
                                           
                                           ) 
                                         
                                       
                                       , 
                                       
                                         l 
                                         
                                           x 
                                           i 
                                         
                                         
                                           ( 
                                           
                                             k 
                                             - 
                                             1 
                                           
                                           ) 
                                         
                                       
                                       , 
                                       
                                         l 
                                         
                                           y 
                                           i 
                                         
                                         
                                           ( 
                                           
                                             k 
                                             - 
                                             1 
                                           
                                           ) 
                                         
                                       
                                     
                                     ) 
                                   
                                 
                               
                             
                             ] 
                           
                         
                       
                       ] 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
     After convergence at step k, the new values for l x     i     (k) , l y     i     (k) , and η i   (k)  are calculated by 
                     [             (     l     x   i       (   k   )       )     2         0           0           (     l     y   i       (   k   )       )     2           ]     =           (       Ψ   ^     i     (   k   )       )     T     ⁢     Σ   i     (   k   )       ⁢     Ψ   i     (   k   )         -     2   ⁢     B   ^     ⁢           ⁢   and               (   11   )                 η   i     (   k   )       =         α   i     (   k   )       ·       (       l     x   i       (   k   )       ·     l     y   i       (   k   )         )       -   1             Σ   j     ⁢       α   j     (   k   )       ·       (       l     x   j       (   k   )       ·     l     y   j       (   k   )         )       -   1                     (   12   )               
under the constraint Σ j η j =1, and where {circumflex over (B)} is the diagonal matrix of estimated imaging PSF correlation lengths and {circumflex over (Ψ)} i   (k)  is found via singular value decomposition of Σ i   (k) .
 
     Equations (9) and (10) are then updated with new values for f(η i ,l x     i   ,l y     i   ) (k)  and Θ (k) . The parameter set Θ (k+1)  is updated via the EM algorithm in the next iteration. The algorithm iterates between the aforementioned update steps until convergence. 
     Upon convergence, the estimates {circumflex over (l)} x     i   , {circumflex over (l)} y     i   , and {circumflex over (η)} i , i=1, . . . , M, are recovered via Equations (11) and (12). Using Equation (4), the variable {circumflex over (θ)} i  is recovered by 
                         θ   ^     i     =       tan     -   1       ⁡     (           Ψ   ^     i     ⁡     (     1   ,   2     )             Ψ   ^     i     ⁡     (     1   ,   1     )         )         ,           (   13   )               
where {circumflex over (Ψ)} i (m,n) is the element on the m-th row and n-th column of matrix {circumflex over (Ψ)} i .
 
     The update equations for the parameters Θ (k)  are as follows: 
                       〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j     =         α   i     (     k   -   1     )       ⁢       G   i     ⁡     (         [         x   j     _     ,       y   j     _       ]     |     [       μ     X   i       (     k   -   1     )       ⁢     μ     Y   i       (     k   -   1     )         ]       ,     Σ   i     (     k   -   1     )         )               Σ     i   =   1     M     ⁢     α   i     (     k   -   1     )       ⁢       G   i     ⁡     (         [         x   j     _     ,       y   j     _       ]     |     [       μ     X   i       (     k   -   1     )       ⁢     μ     Y   i       (     k   -   1     )         ]       ,     Σ   i     (     k   -   1     )         )         ,               (   14   )                       ⁢         m   j     ⁡     (     Θ     (   k   )       )       =     {               n   j     ,     j   =   1     ,   …   ⁢           ,   r                 n   ·         P   j     ⁡     (     Θ   j     (   k   )       )         P   ⁡     (     Θ     (   k   )       )           ,     j   =     r   +   1       ,   …   ⁢           ,   v           ,                 (   15   )                       ⁢         α   i     (   k   )       =         Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j           Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )             ,             (   16   )                       ⁢           [       μ     X   i       (   k   )       ⁢     μ     Y   i       (   k   )         ]     T     =         Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢           〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j     ⁡     [         x   j     _     ,       y   j     _       ]       T           Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j           ,   and             (   17   )                       ⁢         Σ   i     (   k   )       =         Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j     ⁢     K   T     ⁢   K         Σ     j   =   1     v     ⁢       m   j     ⁡     (     Θ     (   k   )       )       ⁢       〈       τ   i     (   k   )       |     Θ   i     (     k   -   1     )         〉     j           ,             (   18   )               
where K=[[  x j   ,  y j   ]−[μ X     i     (k)  μ Y     i     (k) ]].
 
     Referring to update Equation (15), j=1, . . . , r, r+1, . . . , ν refers to the index of the resolution cells over which the mixture function is defined. Indices 1, . . . , r are cells that contain a valid estimate, while indices r+1, . . . , ν define resolution cells with truncated or missing data. Further, n j  is the bin count of the resolution cell with index j. Also, 
             n   ⁢     =   Δ     ⁢       Σ     j   =   1     r     ⁢       n   j     .             
In addition,
 
                 P   j     ⁡     (   Θ   )         P   ⁡     (   Θ   )             
is the bin probability of the area defined by the resolution cell (x j ,y j ):
 
     
       
         
           
             
               
                 
                   
                     
                       
                         P 
                         j 
                       
                       ⁡ 
                       
                         ( 
                         Θ 
                         ) 
                       
                     
                     = 
                     
                       
                         ∫ 
                         
                           
                             x 
                             j 
                           
                           , 
                         
                       
                       ⁢ 
                       
                         
                           ∫ 
                           
                             y 
                             j 
                           
                         
                         ⁢ 
                         
                           
                             f 
                             ⁡ 
                             
                               ( 
                               
                                 x 
                                 , 
                                 
                                   y 
                                   | 
                                   Θ 
                                 
                               
                               ) 
                             
                           
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             ⅆ 
                             x 
                           
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             ⅆ 
                             y 
                           
                         
                       
                     
                   
                   , 
                   where 
                 
               
               
                 
                   ( 
                   19 
                   ) 
                 
               
             
             
               
                 
                   
                     f 
                     ⁡ 
                     
                       ( 
                       
                         x 
                         , 
                         
                           y 
                           | 
                           Θ 
                         
                       
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         1 
                       
                       M 
                     
                     ⁢ 
                     
                       
                         α 
                         i 
                       
                       ⁢ 
                       
                         
                           G 
                           i 
                         
                         ( 
                         
                           x 
                           , 
                           
                             y 
                             | 
                             
                               [ 
                               
                                 
                                   μ 
                                   
                                     x 
                                     i 
                                   
                                 
                                 ⁢ 
                                 
                                   μ 
                                   
                                     y 
                                     i 
                                   
                                 
                               
                               ] 
                             
                           
                           , 
                           
                             
                               ∑ 
                               i 
                             
                             ⁢ 
                             and 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   20 
                   ) 
                 
               
             
             
               
                 
                   
                     P 
                     ⁡ 
                     
                       ( 
                       Θ 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         j 
                         = 
                         1 
                       
                       r 
                     
                     ⁢ 
                     
                       
                         
                           P 
                           j 
                         
                         ⁡ 
                         
                           ( 
                           Θ 
                           ) 
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
           
         
       
     
     Referring now to  FIG. 1 , there is shown a block diagram of a method  10  for implementing sonar image feature segmentation. At block  100 , a complex beamformed SAS image, S(x,y), is divided into N×O overlapping cells. The amount of overlap is a tradeoff between texture parameter resolution and processing speed. Block  200  estimates the ACF parameters for each cell, as described hereinabove and further described hereinafter. 
     Based on the estimated parameters for a cell, the cell texture parameters are stored at block  300 . If there are more cells, as determined at block  400 , method  10  returns to estimate parameters for the next cell (block  500 ). If all cells have been processed, method  10  classifies or segments the image textures (block  600 ) via pattern classification algorithms known to those of skill in the art. Such image feature or texture segmentation finds further use in such fields as target recognition and environmental characterization. 
     Referring now to  FIG. 2 , there is shown a block diagram corresponding to block  200  of  FIG. 1 . At block  210 , initial estimates for the mean intensity, {circumflex over (μ)} σ , the imaging PSF ACF, {circumflex over (R)} h (X, Y), the intensity ACF, {circumflex over (R)} I (X, Y), the PSF variables {circumflex over (β)} x  and {circumflex over (β)} y , and shape parameter {circumflex over (ν)} are obtained. 
     Referring now to  FIG. 3 , there is shown a block diagram corresponding to block  210  of  FIG. 1 . At block  212 , the mean intensity, {circumflex over (μ)} σ , is estimated using the known expression 
                 μ   ^     I     =       1     N   ⁢           ⁢   O       ⁢       ∑     x   i       ⁢       ∑     y   i       ⁢       S   ⁡     (       x   i     ,     y   i       )       .                 
As is known in the art, the imaging PSF ACF, {circumflex over (R)} h (X, Y), can be estimated (block  214 ) via discrete 2-D Fourier transform:
 
                           R   ^     h     ⁡     (     X   ,   Y     )       =             (       μ   ^     I     )       -   1         N   ⁢           ⁢   O       ·     F     -   1         ⁢     {     F   ⁢     {     S   ⁡     (     x   ,   y     )       }     ⁢     *   •     ⁢     conj   ⁡     (     F   ⁢     {     S   ⁡     (     x   ,   y     )       }       )         }         ,           (   22   )               
where F{•} is the discrete Fourier transform operator, conj(•) is the conjugate operator and             is the element-wise array multiplication operator. Similarly, the intensity ACF, {circumflex over (R)} I (X, Y), is estimated at block  216  via:

     
       
         
           
             
               
                 
                   
                     
                       
                         R 
                         ^ 
                       
                       I 
                     
                     ⁡ 
                     
                       ( 
                       
                         X 
                         , 
                         Y 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         1 
                         
                           N 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           O 
                         
                       
                       · 
                       
                         F 
                         
                           - 
                           1 
                         
                       
                     
                     ⁢ 
                     
                       { 
                       
                         F 
                         ⁢ 
                         
                           
                             { 
                             
                               
                                 
                                    
                                   
                                     S 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         x 
                                         , 
                                         y 
                                       
                                       ) 
                                     
                                   
                                    
                                 
                                 2 
                               
                               ⁢ 
                               
                                 * 
                                 • 
                               
                               ⁢ 
                               
                                 conj 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     F 
                                     ⁢ 
                                     
                                       { 
                                       
                                         
                                            
                                           
                                             S 
                                             ⁡ 
                                             
                                               ( 
                                               
                                                 x 
                                                 , 
                                                 y 
                                               
                                               ) 
                                             
                                           
                                            
                                         
                                         2 
                                       
                                       } 
                                     
                                   
                                   ) 
                                 
                               
                             
                             } 
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   23 
                   ) 
                 
               
             
           
         
       
     
     At block  218 , {circumflex over (β)} x  and {circumflex over (β)} y  are estimated via known least-squares matrix computation using the estimates of {circumflex over (R)} h (X, Y) and {circumflex over (μ)} I  obtained at blocks  214  and  212 , respectively. The estimate of shape parameter {circumflex over (ν)} is obtained at block  220  using the known method for estimating parameters of K-distributed clutter. 
     Referring back to  FIG. 2 , parameters are initialized at block  230 , including the parameter set Θ i ={α i , [μ X     i    μ Y     i   ], Σ i } for i=1, . . . , M, where M is the number of Gaussian mixture components. Additionally, index k is set equal to 0, error variables ε 1  and ε 2  are set, and the likelihood function, L(Θ 0 ), is initialized via Equations (10)-(12). As was the case with the overlap amount, the values for ε 1  and ε 2  are a tradeoff between texture parameter resolution and processing speed. 
     At block  240 , index k is incremented and index j is set equal to 1. Block  250  performs calculations for maximizing the likelihood function. Referring now to  FIG. 4 , there is shown a block diagram corresponding to block  250  of  FIG. 2 . Block  252  calculates E j {τ i   (j) ([  x j   ,  y j   ]νΘ i   (j-1) )} via Equation (14), where E j  is the known membership or partition function. Block  254  calculates α i   (j)  via Equation (16). Equation (17) is used to calculate [μ X     i     (j)  μ Y     i     (j) ] at block  256  and Σ i   (j)  is calculated at block  258  via Equation (18). 
     Block  260  determines if error variable ε 2 ≧|L Θ   (j) −L Θ   (j-1) |, where L Θ   (j)  is obtained from Equations 9 and 10 using the values calculated at block  250 . If, not, method  10  increments index j at block  270  and returns to block  250 . If block  260  is yes, {circumflex over (R)} I     N     (k+1) (X, Y) is calculated (block  280 ) via Equation (5) using {circumflex over (R)} I (X, Y), {circumflex over (μ)} σ , and {circumflex over (R)} h (X, Y) from block  210  and η i   (k) , l x     i     (k) , l y     i     (k) , and [μ X     i     (j)  μ Y     i     (j) ] from block  250 . 
     If |{circumflex over (R)} I     N     (k+1) (X, Y)−{circumflex over (R)} I     N     (k) (X, Y)| is greater than error variable δ 1 , as determined at block  290 , method  10  returns to block  240 . (For clarity of illustration, |{circumflex over (R)} I     N     (k+1) (X, Y)−{circumflex over (R)} I     N     (k) (X, Y)| is shown as |{circumflex over (R)} I     N     (k+1) −{circumflex over (R)} I     N     (k) | in  FIG. 2 ). If the difference is less than or equal to ε 1 , method  10  returns the variables {circumflex over (l)} x     i   , {circumflex over (l)} y     i   , {circumflex over (η)} i , {circumflex over (θ)} i , [{circumflex over (μ)} X     i    {circumflex over (μ)} Y     i   ] for i=1, . . . , M, which are stored at block  300  and used in classifying the cell pattern at block  600  of  FIG. 1 . 
     Referring now to  FIGS. 5A-5D , there are depicted SAS image textures corresponding to hardpack sand ( FIG. 5A ), seagrass ( FIG. 5B ), rocky seabottom ( FIG. 5C ), and sand ripple ( FIG. 5D ). Referring also to  FIGS. 6A-6D , there are shown the log intensity ACF estimates associated with the textures depicted in  FIGS. 5A-5D , respectively. Applying method  10  to the estimated texture ACFs in  FIGS. 6A-6D  results in the reconstructed log intensity ACFs illustrated in  FIGS. 7A-7D , respectively. The natural log of the intensity ACFs are displayed in  FIGS. 6A-6D  and  FIGS. 7A-7D  for ease of visual comparison. 
     Referring now to  FIG. 8A , there is illustrated an exemplary image with two synthetic textures having a diagonal boundary between them running from the left upper corner to the right lower corner of  FIG. 8A . The following parameter values are the same for the two synthetic textures: ν=0.0 μ X =μ Y =0, θ=0°, and β x =β y =0.5. The synthetic textures have differing values of l x =l y =3 and l x =l y =9, corresponding to the lower left and upper right of  FIG. 8A , respectively. 
       FIGS. 8B and 8C  show the parameter estimations for the two textures, with estimated {circumflex over (l)} x  shown in  FIG. 8B  and estimated {circumflex over (l)} y  shown in  FIG. 8C . As can be seen in  FIGS. 8B and 8C , there are distinctive differences in the two texture regions, indicating that the ACF parameter values can be used to segment the image into discernable texture regions. 
     What has thus been described is a sonar image feature extractor wherein a parameterized statistical model is fitted to sonar image data to extract image texture features, such as correlation length, periodicity or ripple, orientation and other statistical descriptors. Such features can be used in automatic target recognition, through-the-sensor environmental characterization and texture segmentation schemes. 
     The sonar image texture segmentation method described herein estimates the parameters of an assumed statistical model from the pixel values in the sonar images. The model expands the known correlated K-distribution model into a Gaussian mixture model to allow more degrees of freedom including modeling periodicity. The model adds an angular component to parameterize rotation of the autocorrelation function about the vertical and horizontal axes. The model further provides a method to estimate texture parameters using an EM method for truncated data. 
     Complex, beamformed SAS images from a variety of sources can be used as input. Additional inputs include a processing cell size (N×O), the amount of cell overlap, the number of ACF mixture components, M, and image error thresholds, ε 1  and ε 2 , for the ACF estimate and the EM algorithm for truncated or missing data, respectively. 
     Outputs of the sonar image texture segmentation method include a vector of estimated ACF model parameters that define the SAS image texture. Using the texture feature vectors for segmentation or classification, a labeled texture image is obtained. 
     The parameterized ACF model described herein is linked to realistic assumptions grounded in the physics of the scattering problem. The extracted features have intuitive meaning such as rotation and correlation length. This is in contrast to Markov random field parameters, which have been found to be difficult to understand when taken alone. Further, Markov parameters typically may rely on restrictive simplifying statistical assumptions, such as Gaussianity in some models. Thus, Markov parameters may not be tractably applied to image modeling problems. Additionally, the method described herein can result in significant man-hour savings in analyzing, documenting, archiving and retrieving seabed image data for hydrographic reconnaissance tasks. 
     Obviously many modifications and variations of the present invention may become apparent in light of the above teachings. For example, the estimated autocorrelation function of the sonar image data could be described by measures other than those described herein. For example, measures such as rotation, major axis length, minor axis length, and peak may be used. These measures would be similar to the ones extracted via the invention. However, in relying on a derivation of a valid statistical model from first principles of scattering within resolution cells, the sonar image feature extractor described herein uses a direct link to a statistical distribution. 
     It will be understood that many additional changes in details, materials, and arrangements of parts which have been described herein and illustrated in order to explain the nature of the invention, may be made by those skilled in the art within the principle and scope of the invention as expressed in the appended claims.