Image enhancement method with simultaneous noise reduction, non-uniformity equalization, and contrast enhancement

[0002] The present invention relates generally to the field of digital image processing, and more particularly to a technique for enhancing the appearance of digital images, such as medical diagnostic images, to more clearly render certain image features while maintaining overall image intensity levels.

[0003] A variety of discrete pixel imaging techniques are known and are presently in use. In general, such techniques rely on the collection or acquisition of data representative of each discrete pixel composing an image matrix. In the medical imaging field, for example, several modalities are available for producing the data represented by the pixels, including magnetic resonance techniques, X-ray techniques, and so forth. Depending upon the particular modality employed, the pixel data is detected and encoded, such as in the form of digital values. The values are linked to particular relative locations of the pixels in the reconstructed image.

[0004] The utility of a processed image is often largely dependent upon the degree to which it can be interpreted by users. For example, in the field of medical diagnostics imaging, MRI, X-ray, and other images are most useful when they can be easily understood and compared by an attending physician or radiologist. Moreover, while a number of image processing parameters may control the final image presentation, it is often difficult to determine which of these parameters, or in what combination the parameters may be adjusted to provide the optimal image presentation. Often, the image processing techniques must be adjusted in accordance with empirical feedback from the physician.

[0005] The facility with which a reconstructed discrete pixel image may be interpreted by an observer may rely upon intuitive factors of which the observer may not be consciously aware. For example, a physician or radiologist may seek specific structures or specific features in an image. In medical imaging, such features might include bone, soft tissue or fluids. Such structures may be physically defined in the image by contiguous edges, contrast, texture, and so forth. The presentation of such features often depends heavily upon the particular image processing technique employed for converting the detected values representative of each pixel to modified values used in the final image. The signal processing technique employed can therefore greatly affect the ability of the observer to visualize salient features of interest. Ideally, the technique should carefully maintain recognizable structures of interest, as well as abnormal or unusual structures, while providing adequate textural and contrast information for interpretation of these structures and surrounding background.

[0006] Known signal processing systems for enhancing discrete pixel images suffer from certain drawbacks. For example, such systems may not consistently provide comparable image presentations in which salient features or structures may be easily visualized. Differences in the reconstructed images may result from particularities of individual scanners and circuitry, as well as from variations in the detected parameters (e.g. molecular excitation or received radiation). Differences can also result from the size, composition and position of a subject being scanned. Signal processing techniques employed in known systems are often difficult to reconfigure or adjust, owing to the relative inflexibility of hardware or firmware devices in which they are implemented or to the coding approach employed in software. Finally, known signal processing techniques often employ computational algorithms which are not particularly efficient, resulting in delays in formulation of the reconstituted image or under-utilization of signal processing capabilities.

[0007] Moreover, certain known techniques for image enhancement may offer excellent results for certain systems, but may not be as suitable for others. For example, low, medium and high field MRI systems may require substantially different data processing due to the different nature of the data defining the resulting images. In current techniques completely different image enhancement frameworks are employed in such cases. In addition, current techniques may result in highlighting of small, isolated areas which are not important for interpretation and may be distracting to the viewer. Conversely, in techniques enhancing images by feature structure recognition, breaks or discontinuities may be created between separate structural portions, such as along edges. Such techniques may provide some degree of smoothing or edge enhancement, but may not provide satisfactory retention of textures at ends of edges or lines.

[0008] With the advent of softcopy displays, problems have presented themselves which did not occur, or were less pronounced with conventional technologies. For example, in x-ray images rendered on photographic film, a light box was typically used by radiologists to read the resulting images, with excellent contrast and detail being provided by the photographic film. Where softcopy is desired, however, the same features may be less visible due to the limited dynamic range of the display device, or to the inherent limitations in rendering the detail in a manner that enables features of interest to be distinguished comparably to photographic film. This problem is particularly acute in view of the fact that softcopy displays are increasingly popular due to the archiving of large volumes of images in institutions, and to the facility of accessing the images for display at workstations where the radiologist will page through a large number of examinations during a reading session.

[0009] There is a need, therefore, for an improved technique for enhancing discrete pixel images which addresses these concerns.

[0010] The present invention provides an image enhancement technique designed to respond to these needs. The invention may be applied in a variety of environments, but is particularly well suited to medical diagnostic imaging, where image data is displayed on softcopy displays such as computer or workstation monitors. The technique offers the advantage of being applicable to both new image-viewing systems of stations, as well as to existing systems.

[0011] In a presently preferred embodiment, the technique first performs an image filtering operation, in which the structural elements in the image are ascertained and separated from the non-structural elements. The structural elements are then orientation smoothed to enhance structural features in these regions. The non-structural regions are homogenized to provide a background for the structural features. The structural regions may be sharpened and certain small features that are identified as noise may be smoothed or removed. Original texture can be added to non-structural regions to further aid in the interpretation of structural and non-structural features.

[0012] The technique then performs a non-uniformity equalization operation, in which pixel values are equalized from original values, while maintaining overall appearance of lighter and darker regions corresponding to higher and lower intensity values. Contrast enhancement is then performed based upon equalized values to render detailed structures more discernable.

[0013] The technique is particularly useful in magnetic resonance imaging systems, but may be employed over a wide range of modalities.

[0014] FIG. 1 is a diagrammatical representation of a magnetic resonance imaging system adapted to enhance discrete pixel images of a subject;

[0015] FIG. 2 is a diagram of an exemplary discrete pixel image made up of a matrix of pixels having varying intensities defining structures and non-structures;

[0016] FIG. 3 is a flow chart illustrating a method for processing a discrete pixel image;

[0017] FIG. 4 is a flow chart illustrating a method for identifying structures in an image and enhancing both structural and non-structural features;

[0018] FIG. 5 is a flow chart illustrating a method for identifying structural features in a digital image;

[0019] FIG. 6 is a diagram of elements used in the steps of FIG. 5 for generating gradient components for each discrete pixel of an image;

[0020] FIG. 7 is a gradient histogram of an image used to identify gradient thresholds for dividing structure from non-structure in the image;

[0021] FIG. 8 is a flow chart illustrating a method for selectively eliminating small or noisy regions from the structure definition;

[0022] FIG. 9 is a flow chart illustrating a method for processing structural features identified in the image;

[0023] FIG. 10 is a flow chart illustrating a method of orientation smoothing structure identified in an image;

[0024] FIG. 11 is a flow chart illustrating a method for performing dominant oriented smoothing in the process summarized in FIG. 10;

[0025] FIG. 12 is a diagram of directional indices employed in the orientation smoothing process of FIG. 11;

[0026] FIG. 13 is a flow chart illustrating a method for performing local orientation smoothing through the process of FIG. 10;

[0027] FIG. 14 is a flow chart illustrating homogenization smoothing of non-structural regions of an image;

[0028] FIG. 15 is a flow chart illustrating orientation sharpening of structural regions in an image;

[0029] FIG. 16 is a flow chart illustrating the reintroduction of textural features of non-structural regions in an image.

[0030] FIG. 17 is a flow chart illustrating non-uniformity equalization in the process of FIG. 3;

[0031] FIG. 18 is a diagrammatical representation of matrix transformations performed in the non-uniformity equalization steps of FIG. 17;

[0032] FIG. 19 is a flow chart illustrating contrast enhancement in the process of FIG. 3;

[0033] FIG. 20 is a diagrammatical representation of matrix transformation operations performed in the contrast enhancement steps of FIG. 19; and

[0034] FIG. 21 is a graphical representation of non-linear dynamic range-mapping transformations implemented in the control logic of FIG. 19.

[0035] Referring to FIG. 1, a magnetic resonance imaging 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 subject 16 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.

[0036] 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.

[0037] 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'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.

[0038] It should be noted that, while in the present discussion reference is made to discrete pixel images generated by a magnetic resonance imaging 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, PET systems, and computer 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 coding in application-specific microprocessors, analog circuitry, or a combination of digital and analog circuitry.

[0039] 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 (see FIG. 1) and may vary depending upon such factors as the subject to be the imaged and the resolution desired.

[0040] Referring now to FIGS. 2 and 3, the image 36 is processed according to control logic implemented by signal processing circuit 24 based upon programming code stored within memory circuit 26. The control logic routine, designated generally by reference numeral 60 in FIG. 3, begins at step 54 where image data from signal acquisition circuitry 22 is input into signal processing circuit 24.

[0041] At step 55 in FIG. 3, the image 36 is filtered by separating structural regions 48 from non-structural regions 50. The structural regions 48 are smoothed and sharpened to aid in the interpretation of the structural features, while the non-structural regions 50 are homogenized to provide a background for the interpretation of the image 36. The separation and filtering steps are described in detail below.

[0042] At step 56, a filtered image is corrected for intensity non-uniformities. A preferred method for correcting for intensity non-uniformities is described in detail below. At step 58, a corrected image is further processed by enhancing the local contrast of the image 36. A preferred method for contrast enhancement is described in detail below. Following step 58, the resulting pixel image values are stored in circuit memory 26 (FIG. 1) for eventual reconstruction into image 36 (FIG. 2) for viewing or storage.

[0043] As illustrated in FIG. 4, the image filtering step 55 of FIG. 3 consists of separating structure regions 48 from non-structural regions 50 and performing enhancement operations on the two types of image regions. As illustrated in FIG. 2, image 36 includes structural regions 48, illustrated as consisting of long, contiguous 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 (i.e., number of adjacent pixels), which may be defined as structural regions, or which may be eliminated from the definition of structure in accordance with the techniques described below.

[0044] This process, designated generally by reference numeral 60 in FIG. 3, begins at step 62 with the initialization of parameters employed in the signal enhancement process. This initialization step includes the reading of default and operator-selected values for parameters described in the following discussion, such 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, requiring the operator to select between several parameter choices, such as image matrix size.

[0045] At step 64, signal processing circuit 24 collects and normalizes the raw values acquired for the pixels defining the image 36. In the illustrated embodiment, this step includes reading digital values representative of intensities at each pixel, and scaling these intensity 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. 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 thus processes the raw image data, denoted Ir in FIG. 4, to render a normalized image In 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.

[0046] It should be noted that while reference is made in the present discussion to intensity values within image 36, 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 frequency or color.

[0047] At step 66, signal processing circuit 24 executes a predetermined logic routine for identifying structure 48 within image 36, as defined by data representative of the individual pixels of the image. Exemplary steps for identifying the structure in accordance with the present technique are described below with reference to FIG. 5. At step 68, the structure identified during step 66 is orientation smoothed as summarized below with reference to FIGS. 10-13. 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 In based upon a structure mask Ms defined during step 66. Orientation smoothing carried out in step 68 thus transforms the normalized image to a filtered image If which will be further refined by subsequent processing.

[0048] The filtered image If is further processed as follows. At step 70, signal processing circuit 24 performs homogenization smoothing on non-structural regions of image 36. As described more fully below with reference to FIG. 14, 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. The process of orientation sharpening is described more fully below with reference to FIG. 15. 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 Ifn to provide background reference for the final image. The texture blending process is described below with reference to FIG. 16. Following step 76, the resulting pixel image values are stored in memory circuit 26 for eventual reconstruction into image 36.

[0049] FIG. 5 illustrates steps in control logic for identifying structural regions 48 within image 36 and for eliminating small or noisy isolated regions from the definition of the structural regions. As indicated above, the logic of FIG. 5, summarized as step 66 in FIG. 4, begins with pixel data of the normalized image In.

[0050] At step 79 a blurred or smoothed version of the normalized image is preferably formed. It has been found that by beginning the steps of FIG. 5 with this smoothed image, structural components of the image may be rendered more robust and less susceptible to noise. While any suitable smoothing technique may be employed at step 79, in the present embodiment, a box-car smoothing technique is used, wherein a box-car filter smoothes the image by averaging the value of each pixel with values of neighboring pixels. As will be appreciated by those skilled in the art, computationally efficient method for such filtering may be implemented, such as employing a separable kernel (3 or 5 pixels in length) which is moved horizontally and vertically along the image until each pixel has been processed.

[0051] At step 80, X and Y gradient components for each pixel are computed based upon the smoothed version of the normalized image. While several techniques may be employed for this purpose, in the presently preferred embodiment, 3×3 Sobel modules or operators 102 and 104 illustrated in FIG. 6, are employed. 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.

[0052] With these gradient components thus computed, at step 82 the gradient magnitude Gmag and gradient direction Gdir are computed. In the presently preferred technique, 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 π/2. The values of the gradient magnitudes and gradient directions for each pixel are saved in memory circuit 26.

[0053] Based upon the gradient magnitude values determined at step 82, a gradient histogram is generated as indicated at step 84. FIG. 7 illustrates an exemplary gradient histogram of this type. 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. 7 need not be represented graphically, but may be functionally determined by the signal processing circuitry operating in cooperation with values stored in memory circuitry.

[0054] 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) or the automatic value identified through the process described below may be overridden by the operator to provide specific information in the resulting image.

[0055] As summarized in FIG. 5, the process for identification of the threshold value begins at step 86 by selecting an initial gradient threshold. This initial gradient threshold, designated 116 in FIG. 7 is conveniently set to a value corresponding to a percentile of the global pixel population, such as 30 percent. 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. 7, 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.

[0056] At step 88, a search is performed for edges of the desired structure. The edge search proceeds 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 presently 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.

[0057] At step 92 small or noisy segments identified as potential candidates for structure are iteratively eliminated. Steps in control logic for eliminating these segments are summarized in FIG. 8. Referring to FIG. 8, 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 (see FIG. 5). 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. 8. Each subsequent iteration is preferably performed in an opposite direction (i.e., from top-to-bottom, and from bottom-to-top).

[0058] 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. 8, 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. 5. 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.

[0059] Returning to FIG. 5, 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 presently 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.

[0060] With the desired number of structure pixels thus identified, a final gradient threshold or FGT is determined as illustrated at step 96 in FIG. 5, based upon the histogram 106 as shown in FIG. 7. 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. 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. In the presently preferred embodiment, the FGT value is then scaled by multiplication by a value which may be automatically determined or which may be set by a user. For example, a value of 1.9 may be employed for scaling the FGT, depending upon the image characteristics, the type and features of the structure viewable in the image, and so forth. The use of a scalable threshold value also enables the technique to be adapted easily and quickly to various types of images, such as for MRI data generated in systems with different field strengths, CT data, and so forth.

[0061] Based upon this scaled 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. 8. 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.

[0062] At step 99, certain of the isolated regions may be recuperated to provide continuity of edges and structures. In the present embodiment, for example, if a pixel in the gradient image is above a second gradient threshold, referred to as GFT, and is connected (i.e. immediately adjacent) to a pixel which is above the FGT, the corresponding pixel in the binary image is changed from a 0 value to a value of 1. The value of the GFT may be set to a desired percentage of the FGT, and may be determined empirically to provide the desired degree of edge and structure continuity. This gradient following step is preferably carried out recursively to determine an initial classification of the pixels.

[0063] At step 100 in FIG. 5, 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. Steps in exemplary control logic for implementing the binary rank order filtering of step 100 are illustrated in FIG. 9.

[0064] Referring to FIG. 9, the binary rank order filtering 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 present 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, control proceeds to step 150.

[0065] 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 present 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.

[0066] With the structure of the image thus identified, orientation smoothing of the structure, as indicated at step 68 of FIG. 4, is carried out through logic such as that illustrated diagrammatically in FIG. 10. As shown in FIG. 10, the orientation smoothing of image structure begins with the normalized values for the corresponding structure and may proceed in different manners depending upon the type of smoothing desired. In particular, based upon an operator input designated 160 in FIG. 10, a logical decision block 162 directs signal processing circuit 24 to either dominant orientation smoothing as indicated at reference numeral 164 or local orientation smoothing as indicated at 166. If dominant orientation smoothing is selected, the intensity values for the structural pixels are processed as summarized below with respect to FIG. 11, to generate a binary mask M′. Following iterations of the procedure outlined below with reference to FIG. 11, the values of mask M′ are evaluated at step 168, and smoothing is performed on the structure intensity values by use of multipliers α and β resulting in values which are then summed as indicated at blocks 170, 172 and 174 of FIG. 10 and as summarized in greater detail below.

[0067] To explain the dominant orientation smoothing step of 164, reference is now made to FIG. 11. As illustrated in FIG. 11, the dominant orientation smoothing begins with assigning directional indices to each pixel identified as a structural pixel in the structure mask Ms. In the present embodiment, one of four directional indices is assigned to each structural pixel in accordance with the statistical variances for each pixel, as shown in FIG. 12. As illustrated in FIG. 12, within a local neighborhood 182 surrounding each structural pixel, statistical variances for pixel kernels in four directions are computed by reference to the normalized intensity values of the surrounding pixels. The direction of the minimum variance is selected from the four computed values and a corresponding directional index is assigned as indicated by reference numeral 184 in FIG. 12. In the present embodiment these directional indices are assigned as follows: “1” for 45 degrees; “2” for 135 degrees; “3” for 90 degrees; and “4” for 0 degrees. These steps are summarized as 186 and 188 in FIG. 11. At step 190 a local area threshold value is assigned based upon the image matrix size. In the present embodiment, a local area threshold of 6 is used for 256×256 pixel images, a value of 14.25 is used for 512×512 pixel images, and a value of 23 is used for 1024×1024 pixel images.

[0068] At step 192, a binary mask M′ is initialized with zero values for each pixel. At step 194 a dominant orientation is established for each structural pixel by examining the directional indices set in step 188 within a local neighborhood surrounding each structural pixel. In this process, the directional indices found in the local neighborhood are counted and the pixel of interest is assigned the directional index obtaining the greatest count (or the lowest index located in the case of equal counts).

[0069] In the present embodiment, both the dominant direction and its orthogonal direction are considered to make a consistency decision in the dominant orientation smoothing operation. In terms of FIG. 12, these directions are 1 and 2, or 3 and 4. It has been found that considering such factors substantially improves the robustness of the dominant orientation determination in the sense of being consistent with the human visual system (i.e. providing reconstructed images which are intuitively satisfactory for the viewer).

[0070] The consistency decision made at step 194 may be based upon a number of criteria. In the present embodiment, the image is smoothed along the dominant direction (i.e. the direction obtaining the greatest number of counts in the neighborhood) if any one of the following criteria is met: (1) the number of counts of the orientation obtaining the greatest number is greater than a percentage (e.g. 67%) of the total neighborhood counts, and the orthogonal orientation obtains the least counts; (2) the number of counts of the orientation obtaining the maximum counts is greater than a smaller percentage than in criterion (1) (e.g. 44%) of the total neighborhood counts, and the orthogonal direction obtains the minimum number, and the ratio of the counts of the dominant direction and its orthogonal is greater than a specified scalar (e.g. 5); or (3) the ratio of the dominant direction counts to its orthogonal direction counts is greater than a desired scalar multiple (e.g. 10).

[0071] In the present embodiment, the neighborhood size used to identify the direction of dominant orientation in step 194 is different for the series of image matrix dimensions considered. In particular, a 3×3 neighborhood is used for 256×256 pixel images, a 5×5 neighborhood is used for 512×512 pixel images, and a 9×9 neighborhood is used for 1024×1024 pixel images.

[0072] At step 196, the count determined in the searched neighborhood for each pixel is compared to the local area threshold. If the count is found to exceed the local area threshold, signal processing circuit 24 advances to step 198. At that step, the intensity value for each structural pixel is set equal to the average intensity of a 1×3 kernel of pixels in the dominant direction for the pixel of interest. Subsequently, at step 200, the value of a corresponding location in the binary matrix M′ is changed from 0 to 1. If at step 196, the count is found not to exceed the local area threshold for a particular pixel, the intensity value for the pixel of interest is set equal to a weighted average as indicated at step 202. This weighted average is determined by the relationship:

[0073] weighted avg=(1/1+p)(input)+(p/l+p)(smoothed value);

[0074] where the input value is the value for the pixel of interest at the beginning of routine 164, p is a weighting factor between 1 and 200, and the smoothed value is the average intensity of a 1×3 kernel in the dominant direction of the pixel of interest. From either step 200 or 202, circuit 24 returns to step 168 of FIG. 10.

[0075] Referring again to FIG. 10, at step 168, the values of each pixel in the binary mask M′ are evaluated. If the value is found to equal zero, the corresponding intensity value Ii is multiplied by a weighting factor α at step 170. In the present embodiment, factor α is set equal to 0.45. At block 172 the resulting value is summed with the product of the normalized intensity value for the corresponding pixel and a weighting factor β as computed at step 174. In the present embodiment, the factors α and β have a sum equal to unity, resulting in a value of β equal to 0.55.

[0076] If at step 168 the value for a particular pixel is found to equal 1 in the binary mask M′, control advances to decision block 176. Decision block 176 is also reached following the summation performed at block 172 as described above. In the present embodiment, the foregoing dominant orientation smoothing steps are performed over a desired number of iterations to provide sufficient smoothing and bridging between structural regions. At step 176, therefore, circuit 24 determines whether the desired number of iterations have been completed, and if not, returns to step 164 to further smooth the structural regions. In the present embodiment, the operator may select from 1 to 10 such iterations.

[0077] As noted above, the orientation smoothing can proceed through an alternative sequence of steps for local orientation smoothing as noted at block 166 in FIG. 10. FIG. 13 illustrates exemplary steps in control logic for such local orientation smoothing. As with the dominant orientation smoothing, the local orientation smoothing begins with the normalized intensity values for the structural pixels. At step 202, statistical variances for 1×3 pixel kernels about each structural pixel are calculated for each indexed direction (see FIG. 12) as described above for the dominant orientation smoothing process. At step 204, a ratio of the maximum/minimum statistical variances identified for each pixel in step 202 is computed. At step 206 this ratio for each structural pixel is compared to a parameter R, referred to as a relaxation factor for the local orientation filtering. In the present embodiment, the value of R can be set between 1 and 200. If at step 206 the variance ratio is found to exceed R, local orientation filtering is accomplished as indicated at step 208 by setting the intensity value for the structural pixel of interest equal to an average value for the 1×3 pixel kernel in the direction of the minimum variance. If at step 206 the ratio between the maximum and minimum variances for the pixel of interest is found not to exceed R, no local orientation smoothing is performed and circuit 24 advances to a point beyond step 208. From this point, control returns to block 170 of FIG. 10.

[0078] As illustrated in FIG. 10, at block 170 the intensity value for each structural pixel is multiplied by a weighting factor α, and combined at block 172 with the product of the normalized intensity value for the corresponding pixel and a weighting factor β produced at block 174. As summarized above, at step 176, circuit 24 determines whether the desired number of iterations has been completed and, if not, returns to the local orientation smoothing block 166, to repeat the steps of FIG. 13 until the desired number of iterations is complete. Once the desired iterations have been performed, the filtered image If resulting from the orientation smoothing is further filtered by the processes described below.

[0079] As summarized above with reference to FIG. 4, following orientation smoothing of the structure identified within the image, homogenization smoothing of non-structure is performed. The steps in a process for such homogenization smoothing are summarized in FIG. 14. As shown in FIG. 14, the normalized intensity values for non-structural pixels are considered in this process. At step 210, the mean neighborhood intensity value for each non-structural pixel is computed (taking into account the normalized values of structural pixels where these are included in the neighborhood considered). In the present embodiment, step 210 proceeds on the basis of a 3×3 neighborhood surrounding each non-structural pixel. This mean value is assigned to the pixel of interest and control advances to step 212. At step 212 circuit 24 determines whether a desired number of iterations has been completed. If not, control returns to step 210 for further homogenization of the non-structural pixel intensity values. Once the desired number of iterations has been completed the homogenization smoothing routine of FIG. 14 is exited. In the present embodiment, the operator may set the number of homogenization smoothing iterations from a range of 1 to 10.

[0080] Following homogenization smoothing, the filtered image is further processed by orientation sharpening of the identified structure as mentioned above with regard to FIG. 4, and as illustrated in greater detail in FIG. 15. As shown in FIG. 15, in the present embodiment, the sharpening is performed only for pixel values which are above a preset lower limit, as indicted at decision block 213. This limit, which may be set to a multiple of the FGT (e.g. 2×FGT), thus avoids enhancement of structural pixels which should not be sharpened. If a structural pixel has a value above the limit, the orientation sharpening sequence begins at step 214 where Laplacian values for each such structural pixel are computed in the indexed directions shown in FIG. 12 and described above. The Laplacian values may be computed from the formula

[0081] L(k)=2.0*I(k)−I(k−1)−I(k+1);

[0082] where k is the structural pixel of interest, “k−1” is the pixel preceding the pixel of interest in the indexed direction, and “k+1” is the pixel succeeding the pixel of interest in the indexed direction. It should be noted that the Laplacian values computed at step 214 are based upon the filtered intensity values (i.e., smoothed values for structure). At step 216, the maximum of the four Laplacian values for each structural pixel is then saved to form an edge image Me. At step 218, for each structural pixel, the statistical variances and mean values for a 3×1 pixel kernel are computed in the indexed directions shown in FIG. 12, again using the filtered (i.e., homogenized and smoothed) values for each pixel in the relevant neighborhoods. The direction of minimum variance for each structural pixel is then identified from these values, and the mean value in the direction of minimum variance is saved for each pixel as indicated at step 220. At step 222, the mean value in the direction of minimum variance for each structural pixel is multiplied by a factor γ. In the present embodiment, the value of γ may be set within a range varying from 0.00001 to 2.0. In general, the higher the value of γ selected, the greater the overall sharpness of strong edges in the final image. At step 224, the resulting weighted value is added to the initial filtered value for the corresponding structural pixel. The resulting filtered image values are then further processed as described below.

[0083] Following orientation sharpening of the structural features of the image, the entire image is again renormalized as indicated at step 74 in FIG. 4. While various methods may be used for this renormalization, in the present embodiment the global average pixel intensity in the filtered image following step 72 is computed, and a normalization factor is determined based upon the difference between this average value and the average value prior to the filtration steps described above. The new normalized intensity value for each pixel is then determined by multiplying this normalization factor by the filtered pixel intensity, and adding the global minimum intensity value from the original data to the product.

[0084] The resulting renormalized filtered image, denoted Ifn in FIG. 4, is then further processed to blend back into the non-structural regions texture from the original data. This texture blending step is summarized in FIG. 16. In general, the steps of FIG. 16 tend to add more or less original texture depending upon the gradient magnitude of the non-structural pixels. In particular, at step 226, the gradient magnitude for each non-structural pixel of interest is compared to a threshold value T. In the present embodiment, this threshold is set to a value of 300. If the gradient is found not to exceed the threshold, the pixel intensity value is multiplied by a value “a” at step 228. The resulting product is added at step 228 to the product of the raw intensity value for the pixel (prior to the normalization at step 64 of FIG. 4) multiplied by a value equal to “1-a” at step 232. The resulting weighted average is assigned to the pixel.

[0085] If at step 226, the gradient magnitude value for a pixel is found to exceed the threshold value T, the pixel intensity is multiplied by a factor “b”, as noted at step 234. The resulting product is then added at step 236 to the product of the raw intensity for that pixel and a multiplier equal to “1-b” determined at step 238. In the present embodiment, the value of “b” may be set within a range from 0 to 1, with the value of “a” being set equal to 1.5 times the value of “b”. As will be apparent to those skilled in the art, the weighted averaging performed by the steps summarized in FIG. 16 effectively adds texture back to non-structural regions to provide understandable environment for the structure filtered as described above. By performing the comparison at step 226, the process effectively adds less original texture for pixels having low gradient values, and more original texture for pixels having higher gradient values. Where desired, the values of “a” and “b” may be set so as to increase or decrease this function of the process.

[0086] Referring again to FIG. 3, after the image data is filtered in step 55, the resulting filtered image is further modified by performing non-uniformity equalization on the image data, which involves shrinking the image data matrix, smoothing and equalizing the image data, expanding the data matrix, and constructing a corrected image based on both the processed image data and the original filtered image.

[0087] Exemplary steps in control logic for carrying out the non-uniform equalization step 56 of FIG. 3 are summarized in FIG. 17. The non-uniform equalization procedure is further illustrated through the diagrammatical representation of FIG. 18. At step 250, the input data for a beginning matrix g and parameters referred to through the following logical steps are first read by the image data processing circuit.

[0088] At step 252, values of g_avgand a threshold value T are computed. The value g_avgis the average value of the intensities of pixels comprising the matrix g. The threshold value T is set equal to the average value g_avgmultiplied by a coefficient t, set by the program, or by an operator to a value typically greater than 0, and in the present embodiment to a value of 0.1. At step 254, a reduced image is computed, denoted gshr, by converting the matrix g to the matrix gshr through the application of shrink parameters S1 and S2. Referring to FIG. 18, the original value matrix g, designated by reference numeral 270 in FIG. 18, has dimensions m×n, and that step 254 is converted to the matrix gshr, having dimensions m′×n′, as designated by reference numeral 272. As will be appreciated by those skilled in the art, this operation is performed by averaging non-overlapping pixel neighborhoods of dimensions S1×S2, as indicated at reference numeral 274, to form single pixels 276 in the matrix gshr.

[0089] At step 256, the matrix gshr is expanded for later performance of the smoothing operation at step 260. As illustrated diagrammatically in FIG. 18, the expansion of matrix gshr is performed by expanding each dimension of the matrix gshr by an input parameter L, which is the dimension of a separable box-car filter. As shown in FIG. 18, the pixel values adjacent to boundaries of the matrix gshr are mirrored into adjacent regions of the expanded matrix gshr-exp as indicated at reference numerals 278 and 280. Thus, as shown at reference numerals 282, pixels adjacent to the edge of the matrix gshr are inserted into locations along the corresponding dimension L of the additional regions 278 and 280. In the presently preferred embodiment, values of input parameter L may be assigned between 3 and 41, although other values may be employed.

[0090] At step 258 a threshold mask image gshr-exp-thr is formed. This threshold matrix is formed by application of the following rules:

[0091] gshr-exp-thr=C1*g_avg,

[0092] if gshr-exp is≦T; and

[0093] gshr-exp-thr=0,

[0094] if gshr-exp is≦T;

[0095] where C1 is a fixed constant set by either the program or by the operator. In a presently preferred embodiment, C1 has a value of 0.01.

[0096] Returning to FIG. 17, as indicated at step 260, the matrices gshr-exp and gshr-exp-thr are next smoothed by use of a box-car filter having a length L. As illustrated diagrammatically in FIG. 18, these matrices 280 and 276 are therefore converted to matrices 286 and 288, respectively, by progressive movement of the box-car filter 290 through the expanded matrices. It should be noted that the box-car filter will smooth the image by the neighborhood average first along the rows of the image followed by another neighborhood average along the columns of the image. Computationally, this process is performed by the use of two addition operations and two subtraction operations per pixel, irrespective of the magnitude of the dimension L.

[0097] Alternatively, the smoothing can be done by multiplying gshr-exp and gshr-exp-thr by +1 or −1 depending on whether the sum of the pixel index is even or odd. Next, Fast Fourier Transforms (FFTs) are performed on gshr-exp and gshr-exp-thr. Other transforms, which will readily occur to those skilled in the art, may be used. Next, a low-pass filtering operation using pre-computed Gaussian coefficients is performed. Such filter operation, which provides a pass band having the shape of a Gaussian curve of selected variance, is considered to be well-known in the art. In alternate embodiments, other techniques, which will readily occur to those skilled in the art, may be employed. The inverse transforms are then computed, resulting in filtered gshr-exp and filtered gshr-exp-thr. These are then multiplied by +1 or −1, again depending on whether the sum of the pixel index is even or odd. This multiplication reverses the effect of the earlier similar step.

[0098] Returning to FIG. 17, at step 262, a correction function hshr is computed by dividing the low pass filtered matrices 286 and 288 created by the smoothing operation of step 260 by one another. Specifically, the matrix 286 which is the low pass filtered matrix of the gshr-exp matrix is divided, within the unexpanded bounds, by the corresponding values of matrix 288, which is the low pass filtered matrix formed from gshr-exp-thr. In this operation, the following rule is applied within the unexpanded coordinates of the pixels:

[0099] hshr=max (lpf(gshr-exp), N1)/max (lpf (gshr-exp-thr), N1);

[0100] where max(a, b) is a function which returns the maximum of the values a and b, and N1 is a constant internal parameter such that 0 is less than N1 is less than 1. In a presently preferred embodiment, N1 is set to a value of 0.0001. As illustrated in FIG. 18, this operation renders a matrix of dimensions m′×n′, as designated by reference numeral 292.

[0101] As noted at step 264 in FIG. 17, a correction function image h, of dimensions m×n, as indicated at reference numeral 294, is next computed by bilinear interpolation of the matrix hshr. As will be appreciated by those skilled in the art, this bilinear interpolation provides for expansion of each pixel 296 of the equalization function into a series of adjacent pixels in a neighborhood 298.

[0102] From the correction function image, a corrected image feq is computed as indicated at step 266 of FIG. 17. In this operation, a corrected image is first computed in accordance with the relation:

[0103] f(x, y)=(g(x, y)*h(x, y))/(h(x, y)*h(x, y)+N2),

[0104] where N2 is a parameter entered by the user, typically between 100 and 5,000. It should be noted that the resulting image is a function both of the correction function and of the original matrix g. From the resulting values, the equalized matrix, represented by reference numeral 300 in FIG. 18, is computed in accordance with the relation:

[0105] feq(x, y)=f(x, y)*g_avg/f_avg,

[0106] where the value f_avgis the average of the values of f(x, y) computed in the previous step.

[0107] To complete the non-uniform equalization operation, a final equalized image f is selected as indicated at step 268. In the presently preferred embodiment, a series of logical rules or expressions are applied for selection of individual pixel values either from the matrix feq or from the original matrix g, depending upon relationships between them. Specifically, the presently preferred rules are intended to avoid suppression of higher intensity values from the matrix g, while avoiding equalization of low intensity values. Exemplary expressions employed in the presently preferred embodiment may be summarized as follows:

[0108] f(x, y)=g(x, y)

[0109] if feq(x, y)<g(x, y); OR

[0110] f(x, y)=g(x, y)

[0111] if feq(x, y)<1; OR

[0112] f(x, y)=g(x, y)

[0113] if feq(x, y)>g(x, y) AND g(x, y)<T; OR

[0114] f(x, y)=g(x, y)

[0115] if (maximum pixel intensity in g)*C2>g(x, y); OR

[0116] f(x, y)=feq(x, y),

[0117] where C2 is a predetermined parameter, set to a value of 0.01 in the presently preferred embodiment. As will be noted from the foregoing, and from the diagrammatical representation of FIG. 18, the matrix f designed by the reference numeral 302, is formed by selection in accordance with these criteria from the matrix feq and the original data matrix g.

[0118] The non-uniform equalization process of step 56 is followed by contrast enhancement, as indicated at step 58 in FIG. 3. This process, as described more fully below, serves to avoid unnecessarily high contrast boost at very strong edges and very low contrast areas in a computationally efficient fashion. At the intermediate frequencies contrast is boosted. The contrast boosting process, in the present embodiment, entails shrinking the data matrix, smoothing the image data, generating a gain function, performing spatially independent scaling based upon the gain function, expansion of the data matrix by bilinear interpolation, and unsharp masking of the resulting image data. The unsharp masking step provides a mid-frequency boosted image. A non-linear mapping operation is then performed on the image data to preserve the mid-frequency boosted image data, while translating image data near the upper and lower extremities of the dynamic range to the useful dynamic range of the softcopy display. Finally, the contrast enhancement operations may include a final scaling step to produce the final processed image data.

[0119] FIG. 19 summarizes exemplary control logic for carrying out the contrast enhancement step 58 of FIG. 3. Referring to FIG. 19, and to the diagrammatical representation of these steps in FIG. 20, as a first step 304, data of a matrix r is read, as well as parameters referred to below through the contrast enhancement steps. In the presently preferred embodiment, the data of the matrix r, designated by reference numeral 328 in FIG. 20, is identical to the final equalized image data for matrix f illustrated diagrammatically in FIG. 18. From this image, an average value r_avgis computed, along with a threshold value TR, as indicated at step 306. As in the non-uniform equalization steps, the average value r_avgis the average of all of the values of the matrix r, while the value TR is the product of this average value and a fraction tr which is input by the operator or may be set by the program, typically greater than 0, and in the present embodiment is set to a value of 0.1.

[0120] At step 308 of FIG. 19, a reduced image matrix rshr is computed as indicated by reference numeral 330 in FIG. 20. The matrix rshr has dimensions m″×n″ as determined by averaging of pixels in non-overlapping neighborhoods of dimensions SR1×SR2, indicated by reference numeral 332 in FIG. 20. These non-overlapping neighborhood pixels are then converted to single pixels 334 in the rshr matrix.

[0121] At step 310, the reduced image is expanded to form a matrix rshr-exp, as designated by reference numeral 336 in FIG. 20. As described above with respect to the non-uniform equalization process, this step is performed by expanding the dimensions of the matrix rshr by portions of dimension LR, indicated by reference numeral 338 in FIG. 20. In a presently preferred embodiment, the parameter LR is set to a value of between 3 and 41. As explained above, the expansion process provides pixel values in the expanded regions which are mirrored from the original matrix, as indicated at reference numeral 340 in FIG. 20.

[0122] At step 312 of FIG. 19, a masked image rshr-exp-thr is created by application of the following relations:

[0123] rshr-exp-thr=CR1*r_avg,

[0124] if rshr-exp>TR; AND

[0125] rshr-exp-thr=0

[0126] if rshr-exp≦TR

[0127] where CR1 is a fixed constant, set to a value of 0.01 in the preferred embodiment.

[0128] Returning to FIG. 19, as indicated at step 314, the resulting matrices 336 and 342 are next smoothed by a low pass filtering operation. As illustrated in FIG. 20, this operation proceeds by application of a separable box-car filter on the matrices rshr-exp and rshr-exp-thr to produce resulting matrices 344 and 346, respectively. The box-car filter has dimensions LR, and proceeds through the rows and columns as described above in the non-uniform equalization process.

[0129] At step 316, a gain function sshr is computed based upon both of the low pass filtered matrices. In particular, in the presently preferred embodiment, the matrix sshr, designated 350 in FIG. 20, is generated in accordance with the expression:

[0130] sshr=max(lpf-rshr-exp, NR1)/max(lpf-rshr-exp-thr, NR1),

[0131] where max(a, b) returns the maximum values of a and b, and NR1 is a constant internal parameter set between 0 and 1. In a presently preferred embodiment, the value of NR1 is set to a value of 0.0001.

[0132] Returning to FIG. 19, as indicated at step 318, the gain function is next scaled to form a matrix ssca, designated 252 in FIG. 20. This scaling function is preferably performed in accordance with the relation:

[0133] ssca=ssch*avg(lpf-rshr-exp)/avg(sshr),

[0134] where the values avg are the average values of the indicated matrices.

[0135] At step 320 illustrated in FIG. 19, a contrast enhanced image s is next computed by linear interpolation of the scaled gain function. As illustrated in FIG. 20, the matrix s, which has dimensions m×n, and is designated by reference numeral 354 in FIG. 20, is populated with pixels in neighborhoods 358 derived from each pixel 356 of the matrix ssca.

[0136] At step 322, a mid-frequency boosted image p is created by unsharp masking of the matrix s. This step is preferably performed in accordance with the relation:

[0137] p(x, y)=NR2*s(x, y)+(1-NR2)*(r(x, y)−s(x, y)),

[0138] where the value NR2 is a user supplied parameter such that NR2>1.

[0139] From the mid-frequency boosted image p, designated by reference numeral 306 in FIG. 20, a non-linear dynamic range mapping function is applied to map the dynamic range of matrix p on to a desired dynamic range of a softcopy display device, as indicated at reference numeral 324 in FIG. 19. The non-linear mapping function expands local contrast within upper and lower regions of the dynamic range. Based upon the values of matrix p, then, the function is applied, in the form of a lookup table 364 (see FIG. 20) to form a mapped matrix pnlut, designated by reference numeral 362 in FIG. 20.

[0140] In the presently preferred embodiment, the non-linear mapping function applied at step 324 serves to appropriately map upper and lower regions of the dynamic range of the softcopy display device, while providing different mapping for a mid region. An exemplary mapping relation is illustrated in FIG. 21. As shown in FIG. 21, the non-linear mapping function 336 may be expressed graphically as a relationship between input values along axis 368, and output values along axis 370. In a central region 372, between lower and upper bounds 374 and 376 of the central region, values of the input (i.e. of the matrix p) are mapped in accordance with the relationship identified by the segment of the function 378. Below the lower boundary of the mid region, a different mapping relationship 380 is applied to bring the lower intensity pixels within this region into the effective dynamic range of the softcopy display device. Similarly, in the upper region to the right of upper boundary 376 of the mid region, a different mapping function 382 is applied.

[0141] As will be appreciated by those skilled in the art, the use of the non-linear mapping function of the type illustrated in FIG. 21 permits the dynamic range of the image pixels to be conformed to the dynamic range of the softcopy display device. In specific applications, the imaging systems described above may support a wide variety of display devices, typically computer monitors, having different output dynamic ranges. In such cases, the mapping function applied at step 324 is preferably adapted for the dynamic range of each supported display device. In the image enhancement process described above, then, the dynamic range of the display device may be access, input, or otherwise recognized, and the mapping performed at step 324 adjusted to conform to the display device. Other display devices, of course, will receive different mapping functions adapted for their specific dynamic ranges.

[0142] To complete the contrast enhancement operation, a final contrast boosted image pboost, designated by reference numeral 384 in FIG. 20, is generated from matrix pnlut by scaling in accordance with the following relation:

[0143] pboost(x, y)=pnlut(x, y)*r_avg/pnlunt_avg,

[0144] where the “avg” values are the average values of the respective matrices.

[0145] While the invention may be susceptible to various modifications and alternative forms, specific embodiments have been shown by way of example in the drawings and have been described in detail herein. However, it should be understood that the invention is not intended to be limited to the particular forms disclosed. Rather, the invention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention as defined by the following appended claims.