Medical image enhancement method for hardcopy prints

[0002] This invention relates in general to a method for the enhancement of medical images for hardcopy prints. More particularly, it relates to such a method which compensates for the losses of image high frequency contents due to the limited capability of the output device and which mimics the image tone scale on the hardcopy prints of the trans-illuminated film on a standard diagnostic lightbox.

[0003] As picture archiving and communication systems (PACs) have been widely installed and tested in the clinical environment, a virtual radiology department with multi-facilities has been gradually established in which medical images can be accessed and viewed using softcopy displays anywhere across the entire computer network. In such a situation, hardcopy prints of medical images seem obsolete and unnecessary. However, this is not true and hardcopies will always be needed for various purposes. For example, for referral medical report an attached image can greatly facilitate the descriptions.

[0004] Photographic inkjet printers and direct thermal printers are among the best choices for non-diagnostic, referral printing due to their wide availability, low cost and portability. However, when using such a printer to directly print medical images for the referral purpose, there are several major challenges: (1) reduced image local contrast and blurred details due to the limited spatial resolution of the printer, and (2) image tone scale which is different from trans-illuminated film because of either the limited density dynamic range of the printer or the reduced viewing luminance.

[0005] The spatial frequency bandwidth of both inkjet and direct thermal printers is limited due to the halftone algorithms used in the printer technology. As a result, the printed image look quite blurred and local contrast in the image is significantly reduced. Consequently, important diagnostic details, which are often associated with the high frequency contents of the image, are undesirably suppressed. To compensate for the limited spatial resolution of the printer, it is essential to enhance the image details prior to printing.

[0006] When a medical image is printed on a reflection medium such as paper, or a trans-illuminated medium such as transparency, its tone scale is usually different from its film image. The optical density which can be produced on the medium and the image viewing luminance are two dominant factors in causing the tone scale difference. For the same image, when viewing in an indoor ambient light condition, a reflection print always looks darker in the mid-tone region than the film displayed on a lightbox, even if they are printed using the same optical density. The luminance of a lightbox used in the standard diagnostic environment is approximately 2,500 cd/m², and the typical minimum and maximum optical densities on a film are 0.21 and 3.0, respectively. Therefore, the maximum transmitance luminance produced by the film/lightbox is more than 1500 cd/m²and an luminance dynamic range of 620:1 can be achieved. Under such a high luminance level and wide dynamic range, a typical human observer can tell around 800 Just Noticeable Differences (JNDs) of luminance changes in an image, based on Barten's model of the human visual system (HVS) using a standard target (H. Blume, “The act/nema proposal for a grey-scale display function standard,” SPIE Medical Imaging, Vol. 2707, pp. 344-360, 1996).

[0007] On the other hand, the luminance of an indoor environment with diffusive light condition varies dramatically, from about fifty to several hundred cd/m². The typical minimum and maximum densities on a photographic reflection print are around 0.1 and 2.2, respectively, which generate a dynamic range of 120:1 and a maximum reflectance luminance of 120 cd/m²if a luminance of 150 cd/m²is assumed for a surface with 100% reflectance. Again based on Barten's model there are only about 400 JNDs that can be perceived on the reflection print. When using a photographic printer to print an image on a transparency, the produced density range is usually lower than a film. If such an image is displayed on the lightbox, a less number of JNDs can be perceived on the transparency than on the film, which make the image look flat.

[0008] Other factors associated with the photographic printers, such as quantization of the input data, further reduce the JNDs that can be perceived on a hardcopy print. Photographic printers ordinarily take only 8-bit data input and therefore can print 256 greylevels at most. On the other hand, film printers often use 12-bit data input and are able to generate sufficient greylevels for diagnostic purposes.

[0009] It is preferable to display a medical image on a hardcopy print in a consistent way and with the similar perceptual image quality of film. In the clinical field, many x-ray images are evaluated on films for primary diagnostics even if digital x-ray images exist and softcopy displays are available. Also, images from other modalities, such as Computed Tomography (CT) or Magnetic Resonance Imaging (MRI), are often printed and viewed on films. A number of reasons can be given for this including cost, legal liability, as well as the unquestionable fact that most radiologists acquired their reading skills by viewing hardcopy films. Therefore, it is highly recommended to use film images as the gold standard for visualization and generate hardcopy prints that have similar tone scale to films.

[0010] Important diagnostic details are often associated with certain image features such as edges or textures. Enhancement of these image features usually involves in increasing local contrast and sharpening edges. For this purpose, histogram equalization and unsharp masking are two major classes of techniques.

[0011] A histogram equalization technique based on the whole image content is suitable to enhance the global image contrast instead of local contrast, and it sometimes attenuates the contrast of the image regions which have scarce population. On the other hand, a histogram equalization technique based on the local image content is designed to enhance the local contrast. However, it is computationally expensive (R. C. Gonzalez and R. E. Woods, Digital Image Processing, pp. 173-189, Addison-Wesley Publishing Company, 1993). In addition, histogram equalization techniques do not sharpen edges in order to address the problem caused by the limited spatial resolution of a printer.

[0012] An unsharp masking technique, broadly speaking, consists of decomposing an image into a low frequency component and a high frequency component, manipulating either component individually then combining them together. Since edges and texture in an image are usually the high frequency content, edge sharpening and local contrast enhancement can be achieved by boosting the high frequency component.

[0013] Various approaches have been developed for unsharp masking technique to decide and control the low and high frequency components of an image. The high frequency component can be determined either using a linear or a nonlinear high pass filter. Certain quadratic operators are commonly used as nonlinear filters to conduct the high pass operation (G. Ramponi, N. Strobel, S. K. Mitra, and T-H. Yu, “Nonlinear unsharp masking methods for image contrast enhancement,” Journal of Electronic Imaging, Vol. 5, No. 3, pp. 353-366, 1996). On the other hand, since the high frequency component of an image can be computed as the difference between the original image and a low pass filtered version of that image, a median filter, which is nonlinear, can be used for low pass filtering to obtain the high frequency component (R. C. Gonzalez and R. E. Woods, Digital Image Processing, pp. 191-195, Addison-Wesley Publishing Company, 1993). However, the computation is usually expensive for a nonlinear filter with a large kernel. A convolution operator with positive coefficients around its center and negative coefficients in the outer periphery, for example an Laplacian filter, constitutes a linear high pass filter. Due to the fact that a linear filter is easy to implement and fast in execution, it has been widely used in unsharp masking techniques. In practice, a rectangle operator with all the coefficients being equal to unity is the most commonly used linear low pass filter, and the high frequency component of an image is obtained with the subtraction from the original image of the low frequency component.

[0014] For the unsharp masking techniques, many methods have been developed to combine the high frequency and low frequency components of an image. Most of them leave the low frequency component unchanged and manipulate only the high frequency component. For example, the high frequency component is multiplied by a boost factor a, then is added to the original image, as given by Eq. 1,

[0015] P=O+α·H,  (1)

[0016] where P is the processed image, O is the original image and H is the high frequency component, respectively. Murakami disclosed a method to determine the boost factor a based on the observing distance and the spatial frequency content of the image (U.S. Pat. No. 5,204,919, issued 1990, inventor Murakami). The boost factor can also be decided based on the statistics of the pixel values of the images, as proposed by Mohamoodi (U.S. Pat. No. 4,571,635, issued 1984, A. B. Mohmoodi). To improve the computation efficiency for hardware execution, Joyce used an annular mean and an annular standard deviation to compute the boost factor (U.S. Pat. No. 4,941,190, issued 1988, inventor T. H. Joyce). Morishita et al proposed to divide the entire image into a plurality of regional images, and for each regional image the boost factor is determined based on either the pixel value of the original image or the standard deviation of the regional pixel values (U.S. Pat. No. 4,794,531, issued 1987, inventors K. Morishita, et al). Song has patented two algorithms, in which the boost factor is determined adaptively based on the variance of the surrounding pixel values, the noise power spectrum of the image system, and the patterns in the image (U.S. Pat. No. 4,783,840, issued 1987, inventor W. J. Song; and U.S. Pat. No. 5,038,388, issued 1989, inventor W. J. Song). Hishinuma et al computed the boost factor based on the median value of a group of pixels within the unsharp mask of a predetermined size (U.S. Pat. No. 4,747,052, issued 1985, inventors K. Hishinuma et al). Kato et al used a boost factor that varies with either the original or its low frequency component (U.S. Pat. No. 4,315,318, issued 1979, inventors H. Kato et al; and U.S. Pat. No. 4,317,179, issued 1979, inventors H. Kato et al). There are a handful of other algorithms that further combine the boost factor with the high frequency component together, which is broadly described by Eq. 2,

[0017] P=O+f(H),  (2)

[0018] where f(H) is a function of H. Haraki et al patented an algorithm in which f(H) varies in a certain way with the absolute value of H (U.S. Pat. No. 5,369,572, issued 1992, inventors Haraki et al). Ishida et. al recommended f(H) increase monotonically with H (U.S. Pat. No. 4,387,428, issued 1980, inventors Ishida et al; and U.S. Pat. No. 4,346,409, issued 1980, inventors Ishida et al). All the above mentioned methods (Eqs. 1-2) are effective in sharpening edges and enhancing local contrast, but none of them tried to address the problems caused by the reduced viewing luminance and the limited dynamic range of the output device, i.e., to propose a solution to compress the dynamic range of the image and adjust the image tone scale.

[0019] To compress the dynamic range of an image, the low frequency component can also be processed, which is mathematically described by Eq. 3,

[0020] P=O−β·L,  (3)

[0021] where β is an attenuation factor. Nakazawa et al used a positive β which is less than unity to reduce the dynamic range of the low frequency component, while leaving the high frequency unchanged (U.S. Pat. No. 5,319,719, issued 1992, inventors Nakazawa et al). Further improvement of this method includes combining β and L together, as given by Eq. 4,

[0022] P=O+g(L),  (4)

[0023] where g(L) is a function of L. Nakazawa et al tried a method in which g(L) is a function that decreases as L increases (U.S. Pat. No. 5,454,044, issued 1994, inventor Nakajima). More specifically, Nakajima recommended g(L) should decrease monotonically as the value of L increases (U.S. Pat. No. 5,454,044, issued 1994, inventor Nakajima; and U.S. Pat. No. 5,608,813, issued 1995, inventor Nakajima). Another algorithm by Nakazawa et al uses a very broad function description (Eq. 5),

[0024] P=A(O)+B(L),  (5)

[0025] where A(O) is a function of O only, and B(L) is a function of L only (U.S. Pat. No. 5,604,780, issued 1995, inventors Nakazawa et al). All the methods in Eqs. 3-5 can be used to effectively compress the image dynamic range, however, they do not sharpen edges and enhance local contrast in order to compensate for the limited spatial resolution of the output devices.

[0026] Several early publications also revealed algorithms that process both the low frequency and high frequency components. Tsuchino et al used an algorithm based on Eq. 6 (U.S. Pat. No. 5,493,622, issued 1994, inventors Tsuchino et al),

[0027] P=O+f₁(L₁)+g(O−L₂).  (6)

[0028] This is a good approach in that it uses f₁(L₁) to compress the dynamic range in a low frequency region L₁of the original image, at the same time uses g(O−L₂) to boost a different high frequency region H₂=O−L₂. Muka et al tested their algorithm in which the attenuation of the low frequency and the amplification of the high frequency are based on the input/output luminance dynamic ranges, the spatial frequency bandwidth of the output device and knowledge of the human visual system (U.S. Pat. No. 5,774,599, issued Jun. 30, 1998, inventors Muka et al). However, no adaptive dependency of attenuation/amplification on individual image pixel value is used by this algorithm. Nevertheless, the image tone scale is still one of the problems to be addressed.

[0029] Other prior art on image enhancement and tone scale adjustment includes the one from Varies (F.P.P.D. Varies, “Automatic, adaptive, brightness independent contrast enhancement,” Signal Processing, No. 21, pp. 169-182, 1990). In his work the author used logarithmic and exponential calculations for dynamic range transform and applied adaptive local difference operator in the logarithmic domain for local contrast enhancement. While good image quality is expected, the computation is complicated and may be formidable for a large image of 2500×2048 pixels.

[0030] In view of the above mentioned drawbacks of the related prior art, the objective of the present invention is to provide a computationally efficient method of digital image processing, especially edge sharpening and tone scale adjustment of medical images, for hardcopy prints which (1) are produced by output devices with limited capabilities, such as a limited spatial resolution and dynamic range, and/or (2) are viewed in a limited viewing conditions, such an indoor ambient light environment.

[0031] According to a feature of the present invention, there is provided a method for edge sharpening and tone scale adjustment of digital images compromising the steps of: (1) obtaining a medical image in digital format; (2) conducting a plurality of operations for Human Visual System (HVS)-based edge enhancement with each operation having a uniquely specified frequency amplification function and range, such that edges are more enhanced in the shadow region than in the highlight region; (3) adjusting the image tone scale to match that of a trans-illuminated film: for a reflection print, the local contrast in the shadow region is enhanced while that in the highlight region is suppressed, and for a transparency print, the local contrast in the shadow region is suppressed while that in the highlight region is enhanced; (4) sending the processed image to the output device for hardcopy prints.

[0032] The present invention has the following advantages:

[0033] 1. Edges sharpening of medical images is human visual system based; edges are more enhanced in the shadow region than in the highlight region, therefore a uniformly perceived edge enhancement can be achieved across the entire image while the noise can be suppressed to a desirable level.

[0034] 2. Tone scale adjustment of medical images is human visual system based. An image tone scale similar to that of a film can be achieved while the dynamic range of the input image is compressed simultaneously.

[0035] 3. Edge enhancement and image tone scale can be adjusted based on the observer's preference, diagnostic regions of interest, and the characteristics of the output devices.

[0036] FIG. 1 is a flow chart of illustrating the method of the present invention.

[0037] FIG. 2 is a diagrammatic view defining the nine sub-regions for high and low frequency decomposition of an image.

[0038] FIG. 3 is a graphical view showing mapping curves for image tone scale adjustment of reflection hardcopy prints. The HVS curve is obtained with L_o=2,500 cd/m², l_o=150 cd/m², D_min=0.21, D_max=3.0, d_m=0.1, and d_max=2.2.

[0039] FIG. 4 is a graphical view of mapping curves for image tone scale adjustment of transparency hardcopy prints.

[0040] FIG. 5 is a block diagram of a digital image processor for carrying out the method of the present invention.

[0041] In general, the present invention is to provide a computationally efficient method of digital image processing, especially edge sharpening and tone scale adjustment of medical images, for hardcopy prints which (1) are produced by output devices with limited capabilities, such as a limited spatial resolution and dynamic range, and/or (2) are viewed in a limited viewing conditions, such an indoor ambient light environment. Both edge sharpening and tone scale adjustment proposed in the present invention are HVS based. Edges are more enhanced in the shadow region of an image than in the highlight region, therefore a uniformly perceived edge enhancement can be achieved across the entire image while the noise can be suppressed to a desirable level. By using tone scale adjustment, for reflection print local contrast in the shadow region is enhanced while that in the highlight region is suppressed, and for transparency print local contrast in the shadow region is suppressed while that in the highlight region is enhanced, therefore an image tone scale similar to that of a film can be achieved while the dynamic range of the input image is compressed simultaneously. Furthermore, the edge enhancement and tone scale adjustment can be adjusted based on the observer's preference, diagnostic regions of interest and the characteristics of the output devices.

[0042] The present invention consists of four major stages of operations as shown in FIG. 1:

[0043] 1. obtaining a medical image in digital format (box 10);

[0044] 2. conducting a plurality of operations for HVS-based edge enhancement with each operation having a uniquely specified frequency amplification function and range (box 12);

[0045] 3. for reflection print, adjusting image tone scale to increase the contrast in the shadow region and decrease the contrast in the highlight region, and for transparency print, adjusting image tone scale to decrease the contrast in the shadow region and increase the contrast in the highlight region (box 14);

[0046] 4. sending the processed image to the output device for hardcopy prints (box 16).

[0047] The reflection hardcopy prints are usually viewed in an indoor ambient light condition, and the transparency prints are displayed on a standard diagnostic lightbox. The hardcopies are produced by a printer capable of printing images on either reflection or transparency media, for example, a direct thermal printer, a inkjet printer, a dot matrix printer, a laser printer or etc.

[0048] The method of the present invention is described as a series of operations performed on a digital medical image. The digital image can be formed by a digital image acquisition system. The image acquisition system can be, for example, (1) a standard x-ray screen/film combination which produces an x-ray film image which is processed chemically or thermally and the processed film digitized by a scanner/digitizer; (2) a computed radiography (CR) system where a latent x-ray image is formed in a storage phosphor and a corresponding digital image is produced by reading out the storage phosphor by a CR reader; (3) a direct radiography (DR) system where an x-ray image is directly acquired by the flat-panel detector array; (4) a direct digital acquisition system typically consisting of a phosphor based scintillating screen coupled to an imager (CCD, MOS) through a lens or fiber optic system; and (5) a diagnostic scanner (such as MRI, CT, US, PET) which produces either directly a digital image or an electronic image which is digitized.

[0049] The digital image is processed in an image processor, according to the method of the present invention. The image processor can take the form of a digital computer, such as illustrated in FIG. 5. In such case, one or more of the steps of said method can be carried out using software routines. The image processor can also include hardware or firmware for carrying out one or more of said method steps. Thus, the steps of the method of the invention can be carried out using software, firmware, and hardware, either alone or in any preferable combination.

[0050] As shown in FIG. 5, a digital computer 300 includes a memory 310 for storing digital images, application programs, operating system, etc. Memory 310 can include mass memory (such as a hard magnetic disc or CD ROM), and fast memory (such as RAM). Computer 300 also includes input device 312 (such as a keyboard, mouse, touch screen), display 314 (CRT monitor, LCD), central processing unit 316 (microprocessor). Components 310,312,314,316 are connected together by control/data bus 320. Computer 300 can include a transportable storage medium drive 318 for reading from and/or writing to transportable storage media 324, such as a floppy magnetic disk or writeable optical compact disk (CD). Computer 300 can include a computer network interface 326 for reading from and/or writing to the computer network 330, in which case the network interface 326 is connected to the control/data bus 320. Computer 300 can include a directly connected output device 322 (thermal printer, dot matrix printer, laser printer, ink jet printer) and/or a computer network 300 connected output device 328. Computer 300 can include an input device 332 which is connected to the computer network. Computer 300 can also include a network storage device 334 which is connected to the computer network.

[0051] As used in this application, computer readable storage medium can include, specifically, memory 310 and transportable storage medium 324. More generally, computer storage medium may comprise, for example, magnetic storage media, such as magnetic disk (hard drive, floppy disk) or magnetic tape; optical storage media, such as optical disk, optical tape, or machine readable bar code; solid state electronic storage devices, such as random access memory (RAM), read only memory (ROM); or any other physical device or medium which can be employed to store a computer program.

[0052] HVS (Human Visual Response) Based Edge Sharpening

[0053] Certain image contents, such as edges and textures, are usually associated with the high frequency component of the image. With a suitable amplification of the high frequency component, edge sharpening and local contrast enhancement can be achieved. Depending on the diagnostic regions of interest, the imaging modality used to obtain the image, and the characteristics of the output devices, determination of the image high frequency component and its corresponding amplification function needs to be considered accordingly. Moreover, since human eyes are typically more sensitive to details in the highlight region of an image than in the shadow region, small details should be generally enhanced more strongly in the shadow region than those in the highlight region in order to achieve a uniformly perceived edge enhancement.

[0054] To improve the computation efficiency, in a preferred embodiment of the present invention, the low frequency component of an image is obtained by convoluting the image with a uniform rectangle convolution kernel, and the high frequency component of the image is obtained by subtracting the low frequency component from the original image. Unlike directly taking the frequency transform of the image, the rectangle convolution kernel only provides a means to decompose an image into the high frequency and low frequency components at a certain frequency separation point; it does not provide detailed frequency information, which makes frequency adaptive amplification impossible. To address this drawback, the present invention uses a plurality of steps of high frequency enhancement with each step having a uniquely specified frequency separation point and amplification function, as given by Eq. 7,

[0055] P_n+1=L_n+α_n·H_n,  (7)

[0056] where n≧0 is the index for processing steps, L_nand H_nare the low and high frequency components of image P_nat processing step n (P₀is the original image O), and α_nbeing the amplification function, which will be described hereinafter. The implementation as expressed by Eq. 7 (1) takes the original image, decomposes it into the high and low frequency components, amplifies the high frequency component using function α₀, then combines the low frequency component and the amplified high frequency component to form an intermediate image, (2) takes the intermediate image as the input and conducts edge enhancement using the same procedure as before except using a different frequency separation point and a different amplification function.

[0057] Determination of the amplification function α_nis based on the human visual system, diagnostic regions of interest and the characteristics of the output device. A constant α_nwill equally enhance all the frequency signals beyond the frequency separation point. A large convolution kernel is often used with a constant amplification function α_nto slightly boost some edges that are associated with the range from the mid-frequency to the high-frequency of an image. Since a normal human observer can better differentiate local contrast changes in the highlight region than in the shadow region, edges in the shadow region need to be boosted more than in the highlight region in order to achieve uniformly perceived edge enhancement. On the other hand, since noise inherent in the image is also greatly boosted when enhancing the high frequency of the image, the edge enhancement needs to be conducted adaptively so that the perceived noise level can be controlled within tolerance. Based on this, function α_nshould be monotonically increasing as the image varies from highlight regions to shadow regions. In a preferred embodiment of the present invention, α_nis a linear function of the pixel intensity L_n, for example, $\begin{array}{cc}{\alpha }_n={\alpha }_{\min }+\left({\alpha }_{\max }-{\alpha }_{\min }\right)·\frac{L_n-L_{\min }}{L_{\max }-L_{\min }},& \left(8\right)\end{array}$

[0058] where L_maxand L_minare the maximum and minimum values of L_n(low frequency component), respectively, and α_maxand α_minare the specified maximum and minimum values of α(boost factors).

[0059] The images to be processed sometimes need a very large memory buffer in the computer. For example a 16-bit image with 2,500×2,048 pixels needs over ten megabytes of storage per image. The method as in Eq. 7 is a memory space saving implementation. However, in certain circumstances the high frequency amplification functions may not be easily controlled since they can be accumulated from step to step due to the use of intermediate images. Depending on the computer resources available, in another preferred embodiment of present invention, the edge enhancement is based on Eq. 9, $\begin{array}{cc}{P}_{n+1}=L_0+\sum _{i=0}^n\left[{\alpha }_n·H_n\right],& \left(9\right)\end{array}$

[0060] where n ≧0 is the index, L₀is the low frequency component of the original image O, H_nis the high frequency component of image O at the nth frequency separation point, P_n+1is the final image, and α_nis the amplification function. The implementation as expressed by Eq. 9 needs more memory space but is easy to control. It first conducts a series of decomposition of the original image into a set of high and low frequency components with each decomposition having a different frequency separation point, then combines each individually amplified high frequency components to the first low frequency component to form a processed image.

[0061] The selection of the shape of the uniform rectangle convolution kernel is depending on the physical height and width of the pixels in the image. Usually a square kernel is used if the physical size of the pixel is also square, otherwise, a non-square kernel must be used based on Eq. 10, $\begin{array}{cc}\frac{K_x}{K_y}=\frac{P_y}{P_x},& \left(10\right)\end{array}$

[0062] where P_xand P_yare the pixel width and height, respectively, and K_xand K_yare the kernel width and height, respectively.

[0063] The selection of frequency separation point, or equivalently, the convolution kernel size, depends on the diagnostic regions of interest, imaging modality and the characteristics of the output device. For CR images, a small kernel (≦1.5 mm) with a large amplification can effectively enhance very small edges and details, however, the image noise is also enhanced simultaneously. Therefore, unless the image itself contains fine and important diagnostic details and structures, such as CR images of bones and joints, small kernel size is not recommended. On the contrary, for CT or MRI images, since the image pixel size usually is very large (˜0.78 mm), a 3×3 kernel is equivalent to a physical area of 2.3×2.3 mm²and is enough to enhance the high frequency component. Table 1 shows a list of example parameters for high frequency amplification of some CR, CT and MRI images when a Kodak Digital Science 8650 PS printer is used as the output device. One can find that bone/joint and chest CR images need two steps of frequency boosts, but abdomen CR and CT/MRI images may only need one boost.

[0064] Once the size of the rectangle kernel is determined, the image decomposition can be accomplished using the conventional convolution method. However, the preferred embodiment of the present invention to decompose an image into low and high frequency components is computed using the following method. First given an image, I(x,y) defined for a rectangular region having corners at (0,0), and (x_max, y_max), the rectangle-of-interest is divided into nine sub-regions. These are defined as follows (FIG. 2):

[0065] sub-region 1: (0,0) to (K_x/2, K_y/2);

[0066] sub-region 2: (K_x/2+1, 0) to (x_max−K_x/2−1, K_y/2);

[0067] sub-region 3: (x_max−K_x/2,0) to (x_max, K_y/2);

[0068] sub-region 4: (0, K_y/2+1) to (K_x/2, y_max−K_y/2−1);

[0069] sub-region 5: (K_x/2+1, K_y/2+1) to (x_max−K_x/2−1, y_max−K_y/2−1);

[0070] sub-region 6: (x_max−K_x/2, K_y/2+1) to (x_max, y_max−K_y/2−);

[0071] sub-region 7: (0, y_max−K_y/2) to (K_x/2, y_max);

[0072] sub-region 8: (K_x/2+1, y_max−K_y/2) to (x_max−K_x/2−1, y_max);

[0073] sub-region 9: (x_max−K_x/2, y_max−K_y/2) to (x_max, y_max);

[0074] where K_xis the number of columns in the kernel, K_yis the number of rows in the kernel, x_maxand Y_maxare the number of columns and rows in the image. Each of the regions is processed as follows:

[0075] sub-regions 1, 2 and 3: vertical image sum buffers are created for each pixel in the first row of the rectangle-of-interest. The vertical image sum buffer for each column of the image is initialized to the sum of the first K_y/2 image pixels in that column. Prior to the calculations for each row in the sub-regions 1, 2 and 3 the vertical image sum buffers are augmented by the image value for that pixel in row y+K_y/2.

[0076] sub-region 1: the image accumulator is initialize to the sum of the first K_x/2+1 vertical image sum buffers. The pixel count accumulator is initialized to (K_x/2+1)×(K_y/2+1). The output image at the first pixel is the value of the image accumulator divided by the value of the pixel count accumulator. For pixel x in the each line the value of the x+K_x/2 vertical image sum buffer is added to the image accumulator and K_y/2+1+y is added to the pixel count accumulator. The output image is the value of the image accumulator divided by the value of the pixel count accumulator.

[0077] sub-region 2: at each successive pixel the value of the x+K_x/2 vertical image sum buffer is added and the value of the x−K_x/2−1 vertical image sum buffer is subtracted from the image accumulator. The pixel count accumulator is unchanged. The output image is the value of the image accumulator divided by the value of the pixel count accumulator.

[0078] sub-region 3: at each successive pixel the value of the x−K_x/2−1 vertical image sum buffer is subtracted from the image accumulator and K_y/2+1+y is subtracted from the pixel count accumulator. The output image is the value of the image accumulator divided by the value of the pixel count accumulator.

[0079] sub-regions 4, 5, and 6: the calculations for these regions are the same as those for regions 1, 2, and 3 respectively with the following exception. Prior to the calculations for each row in the sub-regions 4, 5 and 6 the vertical image sum buffers are augmented by image pixel value in row y+K_y/2 and reduced by the image pixel value in row y−K_y/2−1. The pixel count accumulator is adjusted appropriately in each of the three sub-regions.

[0080] sub-regions 7, 8, and 9: the calculations for these regions are the same as those for regions 1, 2, and 3 respectively with the following exception. Prior to the calculations for each row in the sub-regions 7, 8 and 9 the vertical image sum buffers are reduced by the image pixel value in row y−K_y/2−1. The pixel count accumulator is initialized to (y_max−y+K_y/2+1)×(K_x/2+1). The output image at the first pixel is the value of the image accumulator divided by the value of the pixel count accumulator.

[0081] sub-region 7: the image accumulator is initialize to the sum of the first K_x/2+1 vertical image sum buffers. The pixel count accumulator is initialized to (K_x/2+1)×(K_y/2+1). The output image at the first pixel is the value of the image accumulator divided by the value of the pixel count accumulator. For pixel x in the each line the value of the x+K_x/2+1 vertical image sum buffer is added to the image accumulator and the pixel count accumulator is incremented by y+K_y/2+1−y_maxto reflect the total number of pixels contributing to the convolution.

[0082] sub-region 8: at each successive pixel the value of the x+K_x/2 vertical image sum buffer is added and the value of the x−K_x/2−1 vertical image sum buffer is subtracted from the image accumulator. The pixel count accumulator is unchanged. The output image is the value of the image accumulator divided by the value of the pixel count accumulator.

[0083] sub-region 9: at each successive pixel the value of the x−K_x/2−1 vertical image sum buffer is subtracted from the image accumulator and y+K_y/2+1−y_maxis subtracted from the pixel count accumulator. The output image is the value of the image accumulator divided by the value of the pixel count accumulator.

[0084] Image Tone Scale Adjustment

[0085] In this section, the transform functions for tone scale adjustment are derived. When an image is displayed on a reflection print instead of film/lightbox, its luminance dynamic range is compressed dramatically, as is the perceived luminance. To produce the similar image tone scale, the JNDs which are perceived on the film/lightbox must be linearly mapped to the JNDs on the reflection print. This can be accomplished by using Blume's standard display function (H. Blume, “The acr/nema proposal for a grey-scale display function standard,” SPIE Medical Imaging, Vol. 2707, pp. 344-360, 1996) as shown in Eq. 11, $\begin{array}{cc}\log L_I=\frac{a+c·\ln I+e·{\left(\ln I\right)}^2+g·{\left(\ln I\right)}^3+j·{\left(\ln I\right)}^4}{1+b·\ln I+d·{\left(\ln I\right)}^2+f·{\left(\ln I\right)}^3+h·{\left(\ln I\right)}^4+k·{\left(\ln I\right)}^5},& \left(11\right)\end{array}$

[0086] where I is the index (1 to 1023) of the luminance levels L_Iof the JNDs, with L_Iranging from 0.05 to 4000 cd/m², and a=−1.30119, b=−0.025840, c=0.080243, d=−0.10320, e=0.13647, f=0.028746, g=−0.025468, h=−0.0031979, j=0.0013635 and k=0.00012993.

[0087] Since the input data of a photographic printer is always represented by a certain unique optical density on the media, a transform function or a look-up-table (LUT) can be applied to linearize the printer input/output response based on the input/output characteristics of the printer. Therefore the output optical density of the printer can be linearly represented by its input. Furthermore, since a medical image is most likely displayed on a film, there exist means to transform the medical image to its film print density space. For example, a processed computed radiography image by the Kodak QCW Workstation (L. L. Barski, et al, “New automatic tone scale method for computed radiography,” SPIE Medical Imaging, Vol. 3335, No. 18, 1998), or a CT/MRI image after applying a LUT based on the CIElab model of the human visual system. The objective here is to map an image from its film density space to its reflection print density space based on Blume's standard display function. Once the mapping is accomplished, the image will be printed on the reflection print media, which can create similar tone scale to the film. The density transform function is derived as follows.

[0088] Given a film optical density D, the luminance L produced by this density on a lightbox is expressed as

[0089] L(D)=L₀·10^−D  (12)

[0090] where L₀is the luminance of the lightbox with a typical value of 2,500 cd/m²for example. The JND value JND corresponding to L(D) can be directly obtained from Eq. 11 together with its maximum and minimum values, JND_maxand JND_min. Similarly, given a density d on the reflection print, the luminance l produced by this density on the print can be expressed as

[0091] l(d)=l₀·10^−D  (13)

[0092] where l₀is the luminance from a surface with 100% reflection of the ambient light, and a value of 150 cd/m²can be used for example in an indoor condition. The JND value jnd corresponding to l(d) can be obtained from Eq. 11 as well as its maximum and minimum values, jnd_maxand jnd_min. A linear transform is required between the JNDs of the film/lightbox and the JNDs of the reflection print, which can be expressed as $\begin{array}{cc}\frac{\mathrm{jnd}-{\mathrm{jnd}}_{\min }}{{\mathrm{jnd}}_{\max }-{\mathrm{jnd}}_{\min }}=\frac{\mathrm{JND}-{\mathrm{JND}}_{\min }}{{\mathrm{JND}}_{\max }-{\mathrm{JND}}_{\min }}.& \left(14\right)\end{array}$

[0093] Therefore, a mapping function from any density on the film to a density on the reflection print has been established using Eqs. 12-14, i.e.,

[0094] D→L(D)→JND(L)→jnd(l)→l(d)→d.  (15)

[0095] Given the viewing condition and density characteristics of both the film/lightbox and the reflection print, using the transform as shown in Eq. 15 one can map any density on the film to a density on the reflection print in such a way that the image tone scale on both presentation media are similar. This is referred to as image tone scale adjustment. Such a HVS based mapping curve under certain conditions is shown in FIG. 3. It should be mentioned here that images which are represented in luminance space can be transformed as well to the reflection print density space using the mapping in Eq. 15.

[0096] The transform for HVS-based image tone scale adjustment is obtained with some assumptions: (1) the front surface reflection caused by the diffuse light is minimal since the reflection print can be manually adjusted to a favorable orientation, (2) the low luminance area of the reflection print do not fall out of the Weber region, and (3) et al. In practice, these assumptions are not always true. The standard display function as expressed by Eq. 11 is based on Barten's HVS model, which uses a special test target and signal modulation. Since a medical image contains signals of complicated frequency components and modulations, the image may not exactly follows the HVS curve in FIG. 3. However, one can believe that the curve is still a good first order approximation and that its shape should be correct based on the derivations in this section, i.e., to transform a film density to a reflection print density, higher density needs more contrast enhancement than lower density, or equivalently, local contrast in the shadow region needs more boost than that in the highlight region. Mathematically speaking, for any two film densities, D₁and D₂with D₁<D₂, the derivative of the function d=T(D) used for image tone scale adjustment should satisfy Ineq. 16 $\begin{array}{cc}{{\frac{dT\left(D\right)}{dD}|}_{D=\mathrm{D2}}\ge \frac{dT\left(D\right)}{dD}|}_{D=\mathrm{D1}}\ge 0.& \left(16\right)\end{array}$

[0097] Based on the above assumption, one can generate a series of curves (FIG. 3) which satisfy Ineq.16 for image tone scale adjustment, and the selection of one particular curve should be depending on factors including the observer's preference, the viewing condition and the diagnostic regions of interest. The mapping curves in FIG. 3 are created according Eq. 17 $\begin{array}{cc}\frac{d-d_{\min }}{d_{\max }-d_{\min }}=\frac{D-D_{\min }}{D_{\max }-D_{\min }+k\left(D_{\max }-D\right)},& \left(17\right)\end{array}$

[0098] where d_min, d_max, D_minand D_maxare the minimum and maximum densities of the reflection print and film, respectively, k is a non-negative constant and is used to generate different transform curves. Curve 0-9 are created with k being equal to 0.00, 0.15, 0.27, 0.56, 0.73, 0.92, 1.15, 1.41, 1.72 and so on, respectively. In practice, a value of k=0.56 will produce a fairly pleasing image tone scale. Also one can find in FIG. 3 that the HVS curve is pretty close to curve 1, which shows that Eq. 17 is a good mathematical approximation in generating the transform functions for image tone scale adjustment.

[0099] It should be mentioned here that a CIElab based LUT to transform a CT or MRI image to density space can be any of the transform curves in FIG. 3.

[0100] The transforms as shown in FIG. 3 make the mid-tone region of an image on the reflection hardcopy appear brighter. For transparency prints, since the highest density which can be produced by the printer is usually less than 3.0. When such an image is displayed on the lightbox, the image looks brighter and flatter than film. To compensate for the density limitation and to achieve a comparable image tone scale on the transparency, a set of transforms which can increase the mid-tone density therefore make the image appear darker are shown in FIG. 4. Generally speaking, for any two film densities, D₁and D₂with D₁<D₂, the derivative of the function d=T(D) used for image tone scale adjustment should satisfy $\begin{array}{cc}{{\frac{dT\left(D\right)}{dD}|}_{D=\mathrm{D2}}\le \frac{dT\left(D\right)}{dD}|}_{D=\mathrm{D1}}\ge 0.& \left(18\right)\end{array}$

[0101] In FIG. 4, the curves are created using Eq. 17 with k being set to 0.00, −0.15, −0.27, −0.56, −0.73, −0.92, −1.15, −1.41, −1.72 and so on, respectively.

[0102] It can be readily concluded that the transform functions are capable of accomplishing two tasks simultaneously: (1) generating an image on the hardcopy print with similar tone scale to a film, and (2) dynamic range compression of the input image. Since both operations are accomplished in one step using a single LUT, the computation efficiency has been greatly improved by the present invention.

[0103] The invention has been described in detail with particular reference to certain preferred embodiments thereof, but it will be understood that variations and modifications can be effected within the spirit and scope of the invention.

[0104] 300 computer

[0105] 310 memory

[0106] 312 input device

[0107] 314 display

[0108] 316 central processing unit

[0109] 318 output device

[0110] 320 control/data bus

[0111] 322 output device

[0112] 324 transportable storage media

[0113] 326 network interface

[0114] 328 network output device

[0115] 330 computer network

[0116] 332 network input device

[0117] 334 network storage device