Patent Publication Number: US-6701025-B1

Title: Medical image enhancement using iteration to calculate an optimal non-uniform correction function

Description:
RELATED APPLICATIONS 
     This application is a continuation-in-part of U.S. patent application Ser. No. 09/317,416 filed on May 24, 1999 and entitled “Method and Apparatus for Enhancing and Correcting Digital Images”. 
    
    
     BACKGROUND OF THE INVENTION 
     The field of the invention is the enhancement of digital images acquired using various imaging modalities. More particularly, the invention relates to producing higher-quality digital images using a combination of noise-reduction and non-uniformity correction techniques. 
     Various techniques have been developed for acquiring and processing discrete pixel image data. Discrete pixel images are composed of an array or matrix of pixels having varying properties, such as intensity and color. The data defining each pixel may be acquired in various manners, depending upon the imaging modality employed. Modalities in medical imaging, for example, include magnetic resonance imaging (MRI) techniques, X-ray techniques, and ultrasonic techniques. In general, each pixel is represented by a signal, typically a digitized value representative of a sensed parameter, such as an emission from material excited within each pixel region or radiation received within each pixel region. 
     To facilitate interpretation of the image, the pixel values must be filtered and processed to enhance definition of features of interest to an observer. Ultimately, the processed image is reconstituted for displaying or printing. In many medical applications, an attending physician or radiologist will consult the image for identification of internal features within a subject, where those features are defined by edges, textural regions and contrasted regions. 
     Unless further processing is applied to a digital image, the image is likely to have a poor signal to noise ratio (SNR), resulting in blurred or ambiguous feature edges and non-uniformities in spatial intensity. Structures, textures, contrasts, and other image features may be difficult to visualize and compare both within single images and between a set of images. As a result, attending physicians or radiologists presented with the images may experience difficulties in interpreting the relevant structures. 
     With respect to non-uniformity correction, in many areas of imaging including MRI and computed tomography, acquired images are corrupted by slowly varying multiplicative inhomogeneities or non-uniformities in spatial intensity. Such non-uniformities can hinder visualization of the entire image at a given time, and can also hinder automated image analysis. Such inhomogeneity is a particular concern in MRI, when single or multiple surface coils are used to acquire imaging data. The acquired images generally contain intensity variations resulting from the inhomogeneous sensitivity profiles of the surface coil or coils. In general, tissue next to the surface coil appears much brighter than tissue far from the coil. Therefore, in order to optimally display and film the entire image, the signal variation due to the inhomogeneous sensitivity profile of the surface coil needs to be corrected. 
     Several prior art methods either enhance features or correct for non-uniformities, but not both. For example, existing techniques for enhancing features may require operator intervention in defining salient structures, sometimes requiring processing of raw data several times based on operator adjustments before arriving at an acceptable final image. This iterative process is inefficient and requires a substantial amount of human intervention. Other prior art methods have been developed for enhancing features of the image while suppressing noise. For example, in one known method, pixel data is filtered through progressive low pass filtering steps. The original image data is thus decomposed into a sequence of images having known frequency bands. Gain values are applied to the resulting decomposed images for enhancement of image features, such as edges. Additional filtering, contrast equalization, and gradation steps may be employed for further enhancement of the image. 
     While such techniques provide useful mechanisms for certain types of image enhancement, they are not without drawbacks. For example, gains applied to decomposed images can result in inadvertent enhancement of noise present in the discrete pixel data. Such noise, when enhanced, renders the reconstructed image difficult to interpret, and may produce visual artifacts which reduce the utility of the reconstructed image, such as by rendering features of interest difficult to discern or to distinguish from non-relevant information. 
     Prior art methods such as that disclosed in U.S. Pat. No. 5,943,433 have also been employed for correcting non-uniformities, although not simultaneously with the above-described methods for feature enhancement. Prior art methods for correcting for non-uniformities include various intensity correction algorithms which correct surface coil images by dividing out an estimate of the surface coil&#39;s sensitivity profile. One such method is based on the assumption that distortion arising from use of surface coils generally varies slowly over space. In accordance with that prior art method, a low pass filtering operation is applied to the measured or acquired image signal. For this prior art method to be effective, however, the image signal must not contain sharp intensity transitions. Unfortunately, at least in MRI imaging, an air-lipid interface usually contains sharp intensity transitions which violate the basic assumption that the low frequency content in the scene being imaged is solely due to the inhomogeneity distortion from the surface coil&#39;s sensitivity profile. 
     Accordingly, certain prior art hybrid filtering techniques have been developed. Although these techniques have been effective in accounting for external transitions, they have not been particularly effective in accounting for significant internal transitions (e.g., transitions that occur between the edges of an organ or other tissue structure). 
     As stated before, acquired images are corrupted by slowly varying multiplicative non-uniformities. When such images are corrected using prior-art techniques, substantial noise amplification can occur, which hinders the visualization of salient features. Therefore, it is common to use less correction than optimal to prevent noise amplification. Besides using less correction, the image may be pre-filtered to reduce noise. Such pre-filtering, however, can also remove salient features from the image. Thus, the combination of pre-filtering and non-uniformity correction techniques has not been put into practice because the combination of prior-art methods has resulted in less-than-optimal images. 
     In image processing literature, several techniques are described to separately improve the SNR and non-uniformity in images. Many authors have described enhancing SNR in MRI images by spatial domain filtering. Likewise, several articles describe improving the shading by correcting for the non-uniformity in the images. Usually these two operations are treated as though they are disjointed operations. 
     R. Guillemaud and M. Brady have discussed simultaneous correction for noise and non-uniformity in IEEE Transactions in Medical Imaging, Vol. 16, pp. 238-251 (1997). These authors used the anisotropic diffusion-based technique proposed by G. Gerig et al., IEEE Transactions in Medical Imaging, Vol. 11, pp. 222-232 (1992), for noise reduction for both pre- and post-filtering with non-uniformity correction. They concluded that pre-filtering loses essential details in the non-uniformity corrected images. Therefore, Guillemaud and Brady chose to perform post-filtering of non-uniformity corrected images. This decision indicates that prior art methods obtain visibility of important weak structures at the expense of non-linear noise amplification. 
     What is needed is an automated method and apparatus which improves the visual quality of digital images. Particularly needed is an automated method and apparatus which reduce noise while correcting for non-uniformities in the image, resulting in better-quality images than were possible using prior art techniques. 
     BRIEF SUMMARY OF THE INVENTION 
     The present invention includes a method and apparatus for automatically correcting a digital image, and particularly, correcting for non-uniformities in the image. An optimal non-uniformity function h is calculated in an iterative process in which input parameters are changed and the non-uniformity function h is evaluated. The digital image represented as image function g is then corrected using the optimal non-uniformity function h to produce a corrected image function f. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 illustrates a simplified block diagram of a magnetic resonance imaging system which incorporates the present invention; 
     FIG. 2 illustrates a diagram of an exemplary discrete pixel image made up of a matrix of pixels having varying intensities defining structures and non-structures; 
     FIG. 3 illustrates a flowchart of a method for processing a discrete pixel image; 
     FIG. 4 illustrates a flowchart of a method for identifying structural features in a discrete pixel image; 
     FIG. 5 illustrates a diagram of elements or modules used in the steps of FIG. 4 for generating gradient components for each discrete pixel of the image; 
     FIG. 6 illustrates a gradient histogram of an image used to identify gradient thresholds for dividing structure from non-structure in the image; 
     FIG. 7 illustrates a flowchart of a method for selectively eliminating small or noisy regions from the structure definition; 
     FIG. 8 illustrates a flowchart of a method for processing structural features identified in the image by binary rank order filtering; 
     FIG. 9 illustrates a flowchart of a method for correcting for non-uniformities in an image; and 
     FIG. 10 illustrates a flowchart of a portion of the method of FIG.  9 . 
    
    
     DESCRIPTION OF THE PREFERRED EMBODIMENT 
     The method and apparatus of the present invention provides the ability to produce higher-quality digital images in which noise has been reduced and non-uniformities have been corrected. This is accomplished by using a robust scheme which does not remove noise, but rather treats that noise as textural information. Essentially, the image is segmented before performing appropriate actions. Smoothing is done along the structures while sharpening is done across them. In the non-structure regions containing weak structures and noise, a homogenizing smoothing is performed, and a part of the texture is added back to retain the original spatial characteristics in a mitigated form. 
     The method and apparatus of the present invention substantially reduces unnatural noise amplification by mitigating noise before it is amplified. Intensity non-uniformities in the image are corrected using a preferred technique. The result of the pre-filtering and correction method is a visually pleasing uniform image which is easy to visualize and film. 
     Referring to FIG. 1, there is shown a simplified block diagram of an MRI system which incorporates the present invention. MRI system  10  is illustrated as including a scanner  12  coupled to circuitry for acquiring and processing discrete pixel data. Scanner  12  includes a support structure  14  in which a physical subject  16  (e.g., a human being) may be placed for acquiring images representative of internal features, such as tissues, fluids and so forth. Scanner  12  includes an electromagnet arrangement  18  for producing an electromagnetic field in a manner generally known in the art. Excitation and sensing coils  20  are provided within scanner  12  for exciting gyromagnetic materials within subject  16  and for sensing emissions from the materials. 
     Signals sensed by coils  20  are encoded to provide digital values representative of the excitation signals emitted at specific locations within the subject, and are transmitted to signal acquisition circuitry  22 . Signal acquisition circuitry  22  also provides control signals for configuration and coordination of fields emitted by coils  20  during specific image acquisition sequences. Signal acquisition circuitry  22  transmits the encoded image signals to a signal processing circuit  24 . Signal processing circuit  24  executes pre-established control logic routines stored within a memory circuit  26  to filter and condition the signals received from signal acquisition circuitry  22  to provide digital values representative of each pixel in the acquired image. These values are then stored in memory circuit  26  for subsequent processing and display. 
     Signal processing circuit  24  receives configuration and control commands from an input device  28  via an input interface circuit  30 . Input device  28  will typically include an operator&#39;s station and keyboard for selectively inputting configuration parameters and for commanding specific image acquisition sequences. Signal processing circuit  24  is also coupled to an output device  32  via an output interface circuit  34 . Output device  32  will typically include a monitor or printer for generating reconstituted images based upon the image enhancement processing carried out by circuit  24 . 
     It should be noted that, while in the present discussion reference is made to discrete pixel images generated by an MRI system, the signal processing techniques described herein are not limited to any particular imaging modality. Accordingly, these techniques may also be applied to image data acquired by X-ray systems, ultrasound systems, digital radiography systems, PET systems, and computed tomography systems, among others. It should also be noted that in the embodiment described, signal processing circuit  24 , memory circuit  26 , and input and output interface circuits  30  and  34  are included in a programmed digital computer. However, circuitry for carrying out the techniques described herein may be configured as appropriate using application-specific microprocessors, analog circuitry, or a combination of digital and analog circuitry. 
     FIG. 2 illustrates a diagram of an exemplary discrete pixel image made up of a matrix of pixels having varying intensities defining structures and non-structures. FIG. 2 illustrates an exemplary discrete pixel image  36  produced via system  10 . Image  36  is composed of a matrix of discrete pixels  38  disposed adjacent to one another in a series of rows  40  and columns  42 . These rows and columns of pixels provide a pre-established matrix width  44  and matrix height  46 . Typical matrix dimensions may include 256×256 pixels; 512×512 pixels; 1,024×1,024 pixels, and so forth. The particular image matrix size may be selected via input device  28  (FIG. 1) and may vary depending upon such factors as the subject to be imaged and the resolution desired. 
     As illustrated in FIG. 2, image  36  includes structural regions  48 , illustrated as consisting of long, continuous lines defined by adjacent pixels. Image  36  also includes non-structural regions  50  lying outside of structural regions  48 . Image  36  may also include isolated artifacts  52  of various sizes, which may be defined as structural regions, or which may be eliminated from the definition of structure in accordance with the techniques described below. 
     FIG. 3 illustrates a flowchart of a method for processing a discrete pixel image. In a preferred embodiment, the method of FIG. 3 is performed by signal processing circuit  24  (FIG. 1) based upon appropriate programming code stored within memory circuit  26  (FIG.  1 ). 
     The method  58  of FIG. 3 begins at step  60  when pixel data is acquired using methods generally known in the art. In step  62 , parameters employed in the signal enhancement process are initialized. This initialization step includes the reading of default and operator-selected values for parameters described in the following discussion, such as the size of small regions to be eliminated from structure, a “focus parameter” and so forth. Where desired, certain of these parameters may be prompted via input device  28  (FIG.  1 ), requiring the operator to select between several parameter choices, such as image matrix size. 
     At step  64 , signal processing circuit  24  collects and normalizes the raw values acquired for the pixels defining the image  36  (FIG.  2 ). In the illustrated embodiment, this step includes reading digital values representative of intensities at each pixel, and scaling these intensities values over a desired dynamic range. For example, the maximum and minimum intensity values in the image may be determined and used to develop a scaling factor over the full dynamic range of output device  32  (FIG.  1 ). Moreover, a data offset value may be added to or subtracted from each pixel value to correct for intensity shifts in the acquired data. At step  64 , circuit  24  (FIG. 1) thus processes the raw image data to render a normalized image which includes pixel values filtered to span a desired portion of a dynamic range, such as 12 bits, independent of variations in the acquisition circuitry or subject. 
     It should be noted that while reference is made in the present discussion to intensity values within image  36  (FIG.  2 ), the present technique may be used to process such values or other parameters of image  36  encoded for individual pixels  38 . Such parameters might include, for example, frequency or color. 
     At step  66 , signal processing circuit  24  (FIG. 1) executes a predetermined logic routine for identifying structure  48  (FIG. 2) within image  36 , as defined by data representative of the individual pixels of the image. A preferred method of identifying structure in accordance with step  66  is described in detail in conjunction with FIGS. 4-8. 
     Referring again to FIG. 3, the image is enhancement filtered via steps  68 - 76 . In step  68 , the structure identified during step  66  is orientation smoothed. While various techniques may be employed for this orientation smoothing, in the embodiment described, dominant orientation smoothing may be carried out, which tends to bridge gaps between spans of structure, or local orientation smoothing may be employed to avoid such bridging. Step  68  is performed on the normalized image based upon a structure mask defined during step  66 . Orientation smoothing carried out in step  68  on the structural regions of the digital image, where those structural regions are identifiable based on the structure mask. Orientation smoothing thus transforms the normalized image to a filtered image which will be further refined by subsequent processing. 
     The filtered image is further processed as follows. At step  70 , signal processing circuit  24  (FIG. 1) performs homogenization smoothing on non-structural regions of image  36  (FIG.  2 ). Homogenization smoothing carried out in step  70  on the non-structural regions of the digital image, where those non-structural regions are identifiable based on the structure mask. This homogenization smoothing is intended to blend features of non-structural regions into the environment surrounding the structure identified at step  66 . At step  72 , the structure identified at step  66  is orientation sharpened. At step  74  the filtered image is then renormalized based upon the intensity values after filtering and the original normalized intensity range. Finally, at step  76  texture present in non-structural regions of the image are blended back into the renormalized filtered image to provide background reference for the final image. 
     Following step  76 , the filtered image is corrected for intensity non-uniformities in step  78 . A preferred method of correcting for intensity non-uniformities is described in detail in conjunction with FIGS. 9 and 10. The resulting pixel image values are stored in memory circuit  26  (FIG. 1) for eventual display as image  36  (FIG.  2 ). 
     The preceding flowchart illustrates, at a high level, a preferred embodiment of the method which employs the present invention, where steps  64 - 76  result in substantial noise reduction (and thus an improved SNR), and step  78  results in correction of intensity non-uniformities. Two steps in this method, identifying structure (step  66 ) and correcting for intensity non-uniformities (step  78 ), are discussed in detail below. In particular, step  66  is described in detail in conjunction with FIGS. 4-8, and step  78  is discussed in detail in conjunction with FIGS. 9 and 10. 
     FIG. 4 illustrates a flowchart of a method for identifying structural features (e.g., features  48 , FIG. 2) in a discrete pixel image (e.g., image  36 , FIG.  2 ). As indicated above, the logic of FIG. 4, summarized as step  66  in FIG. 3, begins with pixel data of the normalized image, In. 
     At step  80 , X and Y gradient components for each pixel are computed. While several techniques may be employed for this purpose, in the preferred embodiment, 3×3 Sobel modules or operators  102  and  104  illustrated in FIG. 5, are employed. 
     FIG. 5 illustrates a diagram of elements or modules used in the steps of FIG. 4 for generating gradient components for each discrete pixel of the image. As will be appreciated by those skilled in the art, module  102  is used for identifying the X gradient component, while module  104  is used for identifying the Y gradient component of each pixel. In this process, modules  102  and  104  are superimposed over the individual pixel of interest, with the pixel of interest situated at the central position of the 3×3 module. The intensity values located at the element locations within each module are multiplied by the scalar value contained in the corresponding element, and the resulting values are summed to arrive at the corresponding X and Y gradient components. 
     Referring back to FIG. 4, with these gradient components thus computed, at step  82  the gradient magnitude, Gmag, and gradient direction, Gdir, are computed. In the preferred embodiment, the gradient magnitude for each pixel is equal to the higher of the absolute values of the X and Y gradient components for the respective pixel. The gradient direction is determined by finding the Arctangent of the Y component divided by the X component. For pixels having an X component equal to zero, the gradient direction is assigned a value of p/2. The values of the gradient magnitudes and gradient directions for each pixel are saved in memory (e.g., memory circuit  26 , FIG.  1 ). 
     It should be noted that, in alternate embodiments, different techniques may be employed for identifying the X and Y gradient components and for computing the gradient magnitudes and directions. For example, those skilled in the art will recognize that, in place of the Sobel gradient modules  102  and  104  (FIG.  5 ), other modules such as the Roberts or Prewitt operators may be employed. Moreover, the gradient magnitude may be assigned in other manners, such as a value equal to the sum of the absolute values of the X and Y gradient components. 
     Based upon the gradient magnitude values determined at step  82 , a gradient histogram is generated as indicated at step  84 . FIG. 6 illustrates a gradient histogram of an image used to identify gradient thresholds for dividing structure from non-structure in the image. The histogram, designated by reference numeral  106 , is a bar plot of specific populations of pixels having specific gradient values. These gradient values are indicated by positions along a horizontal axis  108 , while counts of the pixel populations for each value are indicated along a vertical axis  110 , with each count falling at a discrete level  112 . The resulting bar graph forms a step-wise gradient distribution curve  114 . Those skilled in the art will appreciate that in the actual implementation the histogram of FIG. 6 need not be represented graphically, but may be functionally determined by the signal processing circuitry operating in cooperation with values stored in memory circuitry. 
     Histogram  106  is used to identify a gradient threshold value for separating structural components of the image from non-structural components. The threshold value is set at a desired gradient magnitude level. Pixels having gradient magnitudes at or above the threshold value are considered to meet a first criterion for defining structure in the image, while pixels having gradient magnitudes lower than the threshold value are initially considered non-structure. The threshold value used to separate structure from non-structure is preferably set by an automatic processing or “autofocus” routine as defined below. However, it should be noted that the threshold value may also be set by operator intervention (e.g. via input device  28 , FIG. 1) or the automatic value identified through the process described below may be overridden by the operator to provide specific information in the resulting image. 
     Referring again to FIG. 4, the process for identifying the threshold value begins at step  86  by selecting an initial gradient threshold (“IGT”). This initial gradient threshold, designated  116  in FIG. 6, is conveniently set to a value corresponding to a percentile of the global pixel population, such as 30 percent, for example. The location along axis  108  of the IGT value  116  is thus determined by adding pixel population counts from the left-hand edge of histogram  106  of FIG. 6, adjacent to axis  110  and moving toward the right (i.e., ascending in gradient values). Once the desired percentile value is reached, the corresponding gradient magnitude is the value assigned to the IGT. 
     Referring back to FIG. 4, at step  88 , a search is performed for edges of the desired structure. The edge search proceeds, in a preferred embodiment, by locating the pixels having gradient magnitudes greater than the IGT value selected in step  86 , and considering a 5×5 pixel neighborhood surrounding the relevant pixels of interest. Within the 5×5 pixel neighborhood of each pixel of interest, pixels having gradient magnitudes above the IGT and having directions which do not differ from the direction of the pixel of interest by more than a predetermined angle are counted. In the preferred embodiment, an angle of 0.35 radians is used in this comparison step. If the 5×5 neighborhood count is greater than a preset number (3 in the present embodiment), the pixel of interest is identified as a relevant edge pixel. At step  90 , a binary mask image is created, wherein pixels identified as relevant edge pixels in step  88  are assigned a value of 1, while all other pixels are assigned a value equal to zero. 
     At step  92  small or noisy segments identified as potential candidates for structure are iteratively eliminated. FIG. 7 illustrates a flowchart of a method for selectively eliminating small or noisy regions from the structure definition. The process begins at step  120  where a binary image is obtained by assigning a value of 1 to pixels having a gradient magnitude value equal to or greater than a desired value, and a value of zero to all other pixels. This binary image or mask is substantially identical to that produced at step  90  (FIG.  4 ). 
     At step  122 , each pixel having a value of 1 in the binary mask is assigned an index number beginning with the upper-left hand corner of the image and proceeding to the lower right. The index numbers are incremented for each pixel having a value of 1 in the mask. At step  124 , the mask is analyzed row-by-row beginning in the upper left by comparing the index values of pixels within small neighborhoods. For example, when a pixel is identified having an index number, a four-connected comparison is carried out, wherein the index number of the pixel of interest is compared to index numbers, if any, for pixels immediately above, below, to the left, and to the right of the pixel of interest. The index numbers for each of the connected pixels are then changed to the lowest index number in the connected neighborhood. The search, comparison, and reassignment then continues through the entire pixel matrix, resulting in regions of neighboring pixels being assigned common index numbers. In the preferred embodiment, the index number merging step of  124  may be executed several times, as indicated by step  126  in FIG.  7 . Each subsequent iteration is preferably performed in an opposite direction (i.e., from top-to-bottom, and from bottom-to-top). 
     Following the iterations accomplished through subsequent search and merger of index numbers, the index number pixel matrix will contain contiguous regions of pixels having common index numbers. As indicated at step  128  in FIG. 7, a histogram is then generated from this index matrix by counting the number of pixels having each index number appearing in the index matrix. As will be apparent to those skilled in the art, each separate contiguous region of pixels having index numbers will have a unique index number. At step  130 , regions represented by index numbers having populations lower than a desired threshold are eliminated from the definition of structure as determined at step  90  of FIG.  4 . In a presently preferred embodiment, regions having a pixel count lower than 50 pixels are eliminated in step  130 . The number of pixels to be eliminated in this step, however, may be selected as a function of the matrix size and the amount and size of isolated artifacts to be permitted in the definition of structure in the final image. 
     Referring back to FIG. 4, with pixels for small segments eliminated from the binary mask created at step  90 , the number of pixels remaining in the binary mask are counted as indicated at step  94 . While the resulting number may be used to determine a final gradient threshold, it has been found that a convenient method for determining a final gradient threshold for the definition of structure includes the addition of a desired number of pixels to the resulting pixel count. For example, in a preferred embodiment, a value of 4,000 is added to the binary mask count resulting from step  92  to arrive at a desired number of pixels in the image structure definition. This parameter may be set as a default value, or may be modified by an operator. In general, a higher additive value produces a sharper image, while a lower additive value produces a smoother image. This parameter, referred to in the present embodiment as the “focus parameter,” may thus be varied to redefine the classification of pixels into structures and non-structures. 
     With the desired number of structure pixels thus identified, a final-gradient threshold (“FGT”) is determined as illustrated at step  96  in FIG. 4, based upon the histogram  106  (FIG.  6 ). In particular, the population counts for each gradient magnitude value beginning from the right-hand edge of histogram  106  are summed moving to the left as indicated by reference number  132  (FIG.  6 ). Once the desired number of structural pixels is reached (i.e., the number of pixels counted at step  94  plus the focus parameter), the corresponding gradient magnitude value is identified as the final gradient threshold  134 . Based upon this final gradient threshold, a new binary mask is defined by assigning pixels having values equal to or greater than the FGT a value of 1, and all other pixels a value of zero. At step  98 , the resulting binary mask is filtered to eliminate small, isolated segments in a process identical to that described above with respect to step  92  and FIG.  7 . However, at step  98 , rather than a four-connected neighborhood, a eight-connected neighborhood (i.e., including pixels having shared edges and corners bounding the pixel of interest) is considered in the index number merger steps. 
     At step  100  in FIG. 4, the feature edges identified through the previous steps, representative of candidate structures in the image, are binary rank order filtered. While various techniques may be employed for this enhancing identified candidate structures, it has been found that the binary rank order filtering provides satisfactory results in expanding and defining the appropriate width of contiguous features used to define structural elements. 
     FIG. 8 illustrates a flowchart of a method for processing structural features identified in the image by binary rank order filtering. The binary rank order filtering process begins at step  140  with the binary mask generated and refined in the foregoing steps. At step  140 , circuit  24  determines whether each pixel in the binary mask has a value of 1. If the pixel found to have a value of 1 in the mask, a neighborhood count is performed at step  142 . In this neighborhood count, pixels in the binary mask having values of 1 are counted within a 3×3 neighborhood surrounding the structural pixel of interest. This count includes the pixel of interest. At step  144 , circuit  24  determines whether the count from step  142  exceeds a desired count, m. In the preferred embodiment, the value of m used at step  144  is 2. If the count is found to exceed the value m, the value of 1 is reassigned to the pixel of interest, as indicated at step  146 . If, however, the count is found not to exceed the value of m, the pixel of interest is assigned the value of 0 in the mask as indicated at step  148 . Following steps  146  and  148 , or if the pixel is found not to have an original value of 1 in the mask at step  140 , the method proceeds to step  150 . 
     At step  150 , circuit  24  reviews the structure mask to determine whether each pixel of interest has a value of 0. If a pixel is located having a value of 0, circuit  24  advances to step  152  to compute a neighborhood count similar to that described above with respect to step  142 . In particular, a 3×3 neighborhood around the non-structure pixel of interest is examined and a count is determined of pixels in that neighborhood having a mask value of 1. At step  154 , this neighborhood count is compared to a parameter, n. If the count is found to exceed the parameter n, the mask value for the pixel is changed to 1 at step  156 . If the value is found not to exceed n, the mask pixel retains its 0 value as indicated at step  158 . In the preferred embodiment, the value of n used in step  154  is 2. 
     Following step  156  or step  158 , the resulting mask, Ms, contains information identifying structural features of interest and non-structural regions. Specifically, pixels in the mask having a value of 1 are considered to identify structure, while pixels having a value of 0 are considered to indicate non-structure. 
     FIGS. 4-8 illustrated a preferred method for identifying structure, where the structure information can then be used to enhancement filter the image in accordance with steps  68 - 76  of FIG.  3 . Another key step to the method of FIG. 3 is the step of correcting for non-uniformities (step  78 ) in the enhanced image. 
     FIG. 9 illustrates a flowchart of a method for correcting for non-uniformities in an image. The standard equation describing the constitution of image data is: 
     
       
         
           g=h*f+n, 
         
       
     
     where g is the image data, h is the non-uniformity function, f is the shading corrected image, and n is the imaging noise. Essentially, when given g, the method of the preferred embodiment determines h and n. 
     The method begins, in step  200 , by reading in image data, g(x, y, z). In a preferred embodiment, the image data has been pre-processed through steps  64 - 76  of FIG.  3 . In alternate embodiments, only some of steps  64 - 76  may have been carried out on the image data before step  200 . 
     In step  202 , the maximum intensity value, Max[g], and the average intensity value, Avg[g], are computed. These values are used in later steps of the method. 
     The next step  206  computes the value of a threshold parameter T. This threshold T is computed by extracting the salient features from the image and constructing an intensity histogram. The intensity histogram is a plot of pixel intensity verses the number of pixels in the extracted salient features which have the intensity. The threshold is set to: 
     
       
           T =(mode+average)/4 
       
     
     where: mode=intensity at the peak of the histogram; 
     average=average of all intensities in the histogram. 
     Determination of threshold T in this manner results in a more robust average intensity of the salient features. 
     An iterative process is now entered in which the optimal value of the non-uniformity function h is calculated. As will be described in detail below, this function h is employed to correct the image g and it is determined in part by two input parameters, a shrink parameter S and a variance parameter V. The values of these input parameters S and V may be manually set, but in the preferred embodiment these parameters are set to initial values of S=32 and V=5 as indicated at process block  207 . The iterative process now to be described calculates h based on these values of S and V, evaluates h, and adjusts the values of S or V as required. The value of h is recalculated using the adjusted parameter values and the parameters S and V are further adjusted until an optimal value of h is obtained. 
     After the parameters S and V are set, SHRUNK[g] is obtained, in step  208 . In accordance with this operation, the pixel array of g is reduced along each edge by the shrink parameters. For a three-dimensional array, the pixel array would be reduced by shrink parameter S along edges respectively parallel to the x, y, and z axes. For example, if g represents a 256×256 array, and SHRUNK[g] represents a 32×32 pixel array, then S=32. 
     In step  210 , THRESH[SHRUNK[g]] is computed. Essentially, the intensity of respective pixels of SHRUNK[g] are compared with the threshold, T. If the intensity of a particular pixel is less than or equal to T, the pixel is assigned the value of zero. Otherwise, it is assigned a value of A*Avg[g], where A is usefully selected to be 0.01. In other embodiments, other values of A could be used. 
     In step  211 , SHRUNK[g] and THRESH[SHRUNK[g]] are multiplied by A 1 , depending on whether the sum of the pixel index is even or odd. The purpose of this multiplication step is so that radial symmetry can be used in the frequency domain. 
     In step  212 , transforms of SHRUNK[g] and THRESH[SHRUNK[g]] are performed. In a preferred embodiment, a Fast Fourier Transform (FFT) is used, although other transforms, which will readily occur to those of skill in the art, may be employed. For example, a Discrete Cosine Transform could be used. 
     Next, a low pass filtering (LPF) process is performed in step  214 . In the preferred embodiment, respective transform components are multiplied by coefficients predetermined in a Gausian filter operation. Such filter operation, which provides a pass band having the shape of a Gausian curve is well-known in the art. The precise shape of this Gausian curve is determined by the input parameter V which is the variance of the filter curve. 
     The inverse transform is then computed in step  216 , resulting in low-pass filtered images LPF[SHRUNK[g]] and LPF[THRESH[SHRUNK[g]]]. Then, in step  217 , LPF[SHRUNK[g]] and LPF[THRESH[SHRUNK[g]]] are multiplied by Å 1 , again depending on whether the sum of the index is even or odd. This multiplication step reverses the effects of step  211 . 
     Next, a maximizing operation is performed in step  218 , wherein respective pixel intensities of the two filtered functions are compared with a small regularization parameter, 1 2  using the equations: 
     
       
         max(LPF[SHRUNK[ g ]], 1 2 ) and max(LPF[THRESH[SHRUNK[ g ]]], 1 2 ). 
       
     
     Usefully, 1 2 =0.0001, although different values could be used. 
     Essentially, the computed pixel intensity is either kept, if it is greater than 1 2 , or else replaced with the value of 1 2 , if 1 2  is greater. This maximizing operation improves numerical stability in subsequent operations by eliminating division by very small or near-zero numbers. This, in turn, reduces noise amplification. 
     A shrunken form of the distortion function, h, can then be determined, in step  220 , from the maximizing operation as follows: 
      SHRUNK[ h ]=max(LPF[SHRUNK[ g ]], 1 2 )/max(LPF[THRESH[SHRUNK[ g ]]], 1 2 ). 
     SHRUNK[h] is then expanded, in step  222 , to provide the non-uniformity function, h. SHRUNK[h] can be expanded, for example, using linear or other interpolation methods. In the case of a three-dimensional array, trilinear interpolation could be used, for example. In the case of a two-dimensional array, bilinear interpolation could be used. 
     Referring still to FIG. 9, the correction function h is evaluated to determine if the optimal function has been calculated. This determination is described below in detail with reference to FIG.  10 . As indicated at decision block  223 , if the optimal correction function h has not been calculated, the system branches and calculates another correction function h after altering one of the input parameters S or V as indicated at process block  225 . The process of calculating the correction function h is repeated with different values of S and V until the optimal correction function h is determined, or the full range of parameter values are tried and the default values are used. 
     Given the optimal non-uniformity function, h, the corrected function, f, is computed in step  224  from the following relationship, accounting for noise: 
     
       
           f=g*h /( h*h+Ψ   1 ), 
       
     
     where Ψ 1  is a regularization parameter derived from the reciprocal of the SNR. 
     It will be seen that the intensity range of f is reduced from the original intensity range of g as a result of the division shown in the above equation. Accordingly, it is necessary to rescale the function f back to the original intensity range, as illustrated by step  226 . In a preferred embodiment, this is achieved by applying the following relations: 
     i) f uniform (x,y)=f(x,y)*Avg(g)/Avg(f), for all (x,y), where f uniform  is the initial non-uniformity corrected image; 
     ii) Construct a binary mask image mask(x,y) such that 
     if f uniform (X,Y)&lt;g(x,y)&lt;T, mask(x,y)=1 
     else if f uniform (x,y)/20&gt;g(x,y)&lt;T, mask(x,y)=1; 
     else if g(x,y)&lt;Max[g]/100, mask(x,y)=1 
     else mask(x,y)=0 
     (Note that pixels with value 1 are foreground pixels and value 0 are background pixels); 
     iii) Set mask image pixels which are 1 to  0  if they are connected with other 1&#39;s to make the connectivity count under a pre-specified number (e.g., 1500); 
     iv) Set mask image pixels which are 0 to 1 if they are connected with other 0&#39;s to make the connectivity count under a pre-specified number (e.g., 1500); 
     v) Perform a binary dilation operation followed by a binary erosion operation to pen any bridges thinner than the chosen structuring element; 
     vi) Set mask image pixels which are 0 to 1 if they are connected with other 0&#39;s to make the connectivity count under a pre-specified number (e.g., 1500); 
     vii) Merge the corrected and uncorrected data using the following steps: 
     f final (x,y)=g(x,y), if mask(x,y)=1 or f uniform (x,y)&lt;1, or 
     f final (x,y)=f uniform (x,y) otherwise. 
     In this step, the final non-uniformity corrected image, f final , is reconstructed from initial non-uniformity corrected image, f uniform , and the input image to the non-uniformity correction process. 
     Referring particularly to FIG. 10, the method for determining if an optimal non-uniformity function h has been calculated performs a number of tests. First, as indicated at process block  250  a number NP, is calculated using the intensity level at the peak of the structural image histogram. The structural image histogram is a count of the number of pixels in the masked image depicting only structures (i.e. no background) at each possible intensity level. The intensity level having the highest number of pixels in this histogram is chosen and the number of pixels NP in the correction function h having values greater than this maximum intensity level are counted. In other words, the number of pixels NP in the correction function h which have values greater than the most predominant brightness level in the structural part of the acquired image are counted. 
     As indicated at process block  252 , another count (NPT) is also made from the correction function h. In this step, the number of pixels NPT in the correction function h having an intensity greater than the calculated threshold T are counted. As indicated at decision block  254 , the first test is to determine the following: 
     
       
         NP&gt; t   1 *NPT, 
       
     
     where t 1  is a constant less than one which is selected by the system based on the class of image acquired. In the preferred embodiment t 1  is set to 0.2. 
     If NP is greater, the optimal correction function h has been calculated and it is used to correct the acquired image g as described in detail above and indicated generally in FIG. 10 at process block  256 . 
     Referring still to FIG. 10, if the value of NP is not large enough, a test is made at decision block  258  to determine if the shrink parameter S is at a preset maximum level, S max . If not, the value of S is doubled, as indicated at process block  260 , and the system branches back to recalculate the correction function h as described in detail above and shown generally in FIG. 10 by process block  262 . The value of S max  is “128” in the preferred embodiment used with a 256×256 pixel image. 
     When the allowable range of shrink parameters S has been tried without satisfying the test in decision block  254 , the variance parameter V is decreased. First, the current value of variance parameter V is compared at decision block  264  with a stored minimum value V min . In the preferred embodiment V min  is set to “2”. If V is not less than this minimum, the variance parameter V is reduced (i.e. divided by two) as indicated at process block  266 . The correction function h is then recalculated using this new variance parameter value and the cycle repeats. 
     If the optimal correction function h is not calculated after the parameters S and V have reached their limits S max  and V min , a test is made at decision block  268  to determine if these limiting values should be used. This is performed by evaluating the following expression: 
     
       
         NP&gt; t   2 *NPT, 
       
     
     where t 2  is a constant value less than one. In the preferred embodiment it is set to 0.05. If NP is greater, the most recently calculated correction function is used to correct the image function g at process block  256 . Otherwise, the S and V parameters are set to their default values as indicated at process block  270 , and the correction function h is recalculated at process  272  and used to correct the image function g. 
     The optimal correction function h is thus calculated in an iterative process in which certain adjustable parameters are automatically and methodically adjusted. While the procedure allows these parameters to be manually set and adjusted, the system will automatically determine the best settings to provide consistently good results.