U.S. patent application number 14/810295 was filed with the patent office on 2016-01-28 for method and system for processing an image.
This patent application is currently assigned to Ramot at Tel-Aviv University Ltd.. The applicant listed for this patent is Ramot at Tel-Aviv University Ltd., Tel HaShomer Medical Research Infrastructure and Services Ltd.. Invention is credited to Yuval BARKAN, Lonia KANELOVITCH, Arie RUNDSTEIN, Miriam SKLAIR-LEVY, Hedva SPITZER, Itzchak YACOV.
Application Number | 20160027162 14/810295 |
Document ID | / |
Family ID | 43532717 |
Filed Date | 2016-01-28 |
United States Patent
Application |
20160027162 |
Kind Code |
A1 |
SPITZER; Hedva ; et
al. |
January 28, 2016 |
METHOD AND SYSTEM FOR PROCESSING AN IMAGE
Abstract
A method of processing an image having a plurality of
picture-elements is disclosed. The method comprises: preprocessing
the image to selectively enhance contrast of some picture-elements
of the image, and employing a companding procedure for the
preprocessed image, thereby providing a processed image. The
contrast enhancement optionally and preferably features a mapping
function having a plateau region for diminishing image data with
intensity levels below an intensity threshold, thereby providing a
preprocessed image.
Inventors: |
SPITZER; Hedva; (Tel-Aviv,
IL) ; BARKAN; Yuval; (Kfar-Sirkin, IL) ;
KANELOVITCH; Lonia; (Holon, IL) ; YACOV; Itzchak;
(Ramat-Gan, IL) ; RUNDSTEIN; Arie; (Herzlia,
IL) ; SKLAIR-LEVY; Miriam; (Mevasseret Zion,
IL) |
|
Applicant: |
Name |
City |
State |
Country |
Type |
Ramot at Tel-Aviv University Ltd.
Tel HaShomer Medical Research Infrastructure and Services
Ltd. |
Tel-Aviv
Ramat-Gan |
|
IL
IL |
|
|
Assignee: |
Ramot at Tel-Aviv University
Ltd.
Tel-Aviv
IL
Tel HaShomer Medical Research Infrastructure and Services
Ltd.
Ramat-Gan
IL
|
Family ID: |
43532717 |
Appl. No.: |
14/810295 |
Filed: |
July 27, 2015 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
13501110 |
Jun 28, 2012 |
9092878 |
|
|
PCT/IL2010/000831 |
Oct 12, 2010 |
|
|
|
14810295 |
|
|
|
|
61272604 |
Oct 13, 2009 |
|
|
|
Current U.S.
Class: |
382/131 |
Current CPC
Class: |
H04N 1/387 20130101;
G06T 5/008 20130101; H04N 1/407 20130101; G06T 2207/10084 20130101;
G06T 5/009 20130101; G06T 2207/10072 20130101; G06T 2207/30068
20130101; G06T 2207/10124 20130101; G06T 2207/10081 20130101; G06T
5/40 20130101 |
International
Class: |
G06T 5/40 20060101
G06T005/40; G06T 5/00 20060101 G06T005/00 |
Claims
1. A method of processing an image having a plurality of
picture-elements, comprising: removing isolated objects from the
image; preprocessing the image to diminish image data with
intensity levels below an intensity threshold, thereby providing a
preprocessed image; employing a companding procedure for said
preprocessed image, thereby providing a processed image; and
outputting said processed image to a computer readable medium or a
display device.
2. The method of claim 1, wherein said removal of said isolated
objects features a connected components labeling procedure.
3. The method of claim 1, wherein said intensity threshold is an
average intensity over all picture-elements of the image having
intensity levels within a predetermined intensity range.
4. The method of claim 1, wherein the image is an image of a
section of a living body and said intensity threshold is correlated
with intensity levels characterizing fatty tissue present in said
body section.
5. The method of claim 1, wherein said mapping function comprises a
monotonic region for intensity levels above said intensity
threshold.
6. The method of claim 1, wherein said companding procedure
comprises, for each picture element of said preprocessed image: (i)
defining a surrounding region of picture elements and a remote
region of picture elements; (ii) calculating at least one
adaptation weight; and (iii) using said at least one adaptation
weight and intensity levels of each picture element of said
surrounding and said remote regions, for assigning a new intensity
level to said picture element, said new intensity level being in a
dynamic range which is smaller than the dynamic range of said
preprocessed image.
7. The method of claim 6, wherein said adaptation weight describes,
at least in part, contrast induction between said surrounding
region and said remote region.
8. The method of claim 7, wherein said adaptation weight features a
plurality of different resolutions, and wherein said contrast
induction is calculated separately for each resolution.
9. The method of claim 8, wherein for each resolution, said
contrast induction is calculated as a subtraction between intensity
levels integrated over said surrounding region and intensity levels
integrated over a central region.
10. The method of claim 1, further comprising employing anisotropic
diffusion procedure for reducing noise.
11. The method of claim 1, wherein said image is a mammography
image.
12. The method of claim 11, wherein said processed image exhibits
Cooper's ligaments.
13. The method of claim 11, further comprising identifying and
highlighting Cooper's ligaments in said processed image.
14. The method of claim 1, wherein said image is selected from the
group consisting of a CT image, an X-ray image other than a CT
image, a gamma image, a PET image, a SPECT image and a thermal
image.
15. Non-transitory computer-readable medium having stored thereon a
computer program, wherein said computer program comprising code
means that when executed by a data processing system carry out the
method according of claim 1.
16. A system for processing an image having a plurality of
picture-elements, comprising: an electronic input for receiving the
image; a preprocessing unit configured for preprocessing the image
to remove isolated objects from the image, and to diminish image
data with intensity levels below an intensity threshold; and a
companding unit which compands said preprocessed image.
17. The system of claim 16, further comprising an imaging system
for acquiring the image, wherein said input unit receives the image
from said imaging system.
18. The system of claim 16, wherein said preprocessing unit is
configured for removing isolated objects from the image, prior to
said selective contrast enhancement.
19. The system of claim 16, further comprising a noise reduction
unit configured for employing anisotropic diffusion procedure.
20. The system of claim 16, wherein said image is a mammography
image and wherein the system comprises a functionality for
highlighting Cooper's ligaments in said processed image.
Description
RELATED APPLICATIONS
[0001] This application is a continuation of U.S. patent
application Ser. No. 13/501,110 filed on Jun. 28, 2012, which is a
National Phase of PCT Patent Application No. PCT/IL2010/000831
having International filing date of Oct. 12, 2010, which claims the
benefit of priority of U.S. Provisional Patent Application No.
61/272,604 filed on Oct. 13, 2009. The contents of the above
applications are hereby incorporated by reference as if fully set
forth herein.
FIELD AND BACKGROUND OF THE INVENTION
[0002] The present invention, in some embodiments thereof, relates
to image processing and, more particularly, but not exclusively, to
a method and system for companding an image.
[0003] Commercially available imaging devices are known to acquire
image information across a wide dynamic range of several orders of
magnitude. Typically, although at the time of image capture the
acquired dynamic range is rather large, a substantial portion of it
is lost once the image is digitized, printed or displayed. For
example, most images are digitized to 8-bits (256 levels) per
color-band, i.e., a dynamic range of about two orders of magnitude.
The problem is aggravated once the image is transferred to a
display or a print medium which is often limited to about 50 levels
per color-band.
[0004] The motivation for developing imaging devices capable of
capturing high dynamic range images is explained by the enormous
gap between the performances of the presently available devices and
the ability of the human visual system to acquire detailed
information from a high dynamic range scene. Specifically, the
human visual system, which is capable of acquiring a dynamic range
of several orders of magnitude, can easily recognize objects in
natural light having a high dynamic range. On the other hand, high
quality images suffer, once displayed on a screen or printed as a
hard copy, from loss of information. For example, when the
difference in intensities between a shaded object and its
illuminated surroundings reaches a dynamic range of 2 orders of
magnitudes, presently available display devices may not be able to
recognize it.
[0005] High dynamic range images are also employed in the medial
diagnostic field, where digital imaging systems have become
increasingly important and often preferred over conventional
techniques.
[0006] For example, in current digital X-ray imaging systems,
radiation from a source is directed toward a subject, and a portion
of the radiation passes through the subject and impacts a detector.
The surface of the detector converts the radiation to light photons
which are sensed. The detector is divided into a matrix of discrete
picture elements, and encodes output signals based upon the
quantity or intensity of the radiation impacting each pixel region.
Because the radiation intensity is altered as the radiation passes
through the patient, the images reconstructed based upon the output
signals provide a projection of the patient's tissues similar to
those available through conventional photographic film
techniques.
[0007] However, the dynamic range of the digital detector may be
different from that of a display device, and some image data is
lost during the transmission from the detector to the display, and
the image data is oftentimes adapted to that of the display. For
example, in a typical CT scan, the attenuation coefficient
typically ranges from -1000 (air) to 2500 (bones), namely a dynamic
range of more than three orders of magnitude. Additional background
art includes U.S. Published Application No.
[0008] 20040165086, International Patent Publication Nos.
WO2004/075535 and WO2009/081394, the contents of which are hereby
incorporated by reference.
SUMMARY OF THE INVENTION
[0009] According to an aspect of some embodiments of the present
invention there is provided a method of processing an image having
a plurality of picture-elements. The method comprises:
preprocessing the image to selectively enhance contrast of some
picture-elements of the image using a mapping function having a
plateau region for diminishing image data with intensity levels
below an intensity threshold, thereby providing a preprocessed
image. The method further comprises employing a companding
procedure for the preprocessed image, thereby providing a processed
image, and outputting the processed image to a computer readable
medium or a display device.
[0010] According to an aspect of some embodiments of the present
invention there is provided a system for processing an image having
a plurality of picture-elements. The system comprises: an input
unit for receiving the image, and a preprocessing unit configured
for preprocessing the image to selectively enhance contrast of some
picture-elements of the image using a mapping function having a
plateau region for diminishing image data with intensity levels
below an intensity threshold, hence to provide a preprocessed
image. The system further comprises a companding unit which
compands the preprocessed image.
[0011] According to some embodiments of the invention the invention
the system further comprises an imaging system for acquiring the
image, wherein the input unit receives the image from the imaging
system.
[0012] According to an aspect of some embodiments of the present
invention there is provided a computer-readable medium having
stored thereon a computer program. The computer program comprises
code means that when executed by a data processing system, carry
out the processing operation described herein.
[0013] According to some embodiments of the invention the method
further comprises removing isolated objects from the image, prior
to the selective contrast enhancement.
[0014] According to some embodiments of the invention the removal
of the isolated objects features a connected components labeling
procedure.
[0015] According to some embodiments of the invention the intensity
threshold is an average intensity over all picture-elements of the
image having intensity levels within a predetermined intensity
range.
[0016] According to some embodiments of the invention the image is
an image of a section of a living body and the intensity threshold
is correlated with intensity levels characterizing fatty tissue
present in the body section.
[0017] According to some embodiments of the invention the mapping
function comprises a monotonic region for intensity levels above
the intensity threshold.
[0018] According to some embodiments of the invention the
companding procedure comprises, for each picture element of the
preprocessed image: (i) defining a surrounding region of picture
elements and a remote region of picture elements; (ii) calculating
at least one adaptation weight; and (iii) using the at least one
adaptation weight and intensity levels of each picture element of
the surrounding and the remote regions, for assigning a new
intensity level to the picture element, the new intensity level
being in a dynamic range which is smaller than the dynamic range of
the preprocessed image.
[0019] According to some embodiments of the invention the
adaptation weight describes, at least in part, contrast induction
between the surrounding region and the remote region.
[0020] According to some embodiments of the invention the
adaptation weight features a plurality of different resolutions,
and wherein the contrast induction is calculated separately for
each resolution.
[0021] According to some embodiments of the invention for each
resolution, the contrast induction is calculated as a subtraction
between intensity levels integrated over the surrounding region and
intensity levels integrated over a central region.
[0022] According to some embodiments of the invention the method
comprises employing anisotropic diffusion procedure for reducing
noise. According to some embodiments of the invention the system
comprises a noise reduction unit configured for employing
anisotropic diffusion procedure.
[0023] According to some embodiments of the invention the image is
a mammography image. According to some embodiments of the invention
the processed image exhibits Cooper's ligaments. According to some
embodiments of the invention the method comprises identifying and
highlighting Cooper's ligaments in the processed image. According
to some embodiments of the invention the system comprises a
functionality for highlighting Cooper's ligaments in the processed
image.
[0024] According to some embodiments of the invention the image is
an X-ray image other than a CT image.
[0025] According to some embodiments of the invention the image is
a gamma image.
[0026] According to some embodiments of the invention the image is
a PET image.
[0027] According to some embodiments of the invention the image is
a SPECT image.
[0028] According to some embodiments of the invention the image is
a thermal image.
[0029] Unless otherwise defined, all technical and/or scientific
terms used herein have the same meaning as commonly understood by
one of ordinary skill in the art to which the invention pertains.
Although methods and materials similar or equivalent to those
described herein can be used in the practice or testing of
embodiments of the invention, exemplary methods and/or materials
are described below. In case of conflict, the patent specification,
including definitions, will control. In addition, the materials,
methods, and examples are illustrative only and are not intended to
be necessarily limiting.
[0030] Implementation of the method and/or system of embodiments of
the invention can involve performing or completing selected tasks
manually, automatically, or a combination thereof. Moreover,
according to actual instrumentation and equipment of embodiments of
the method and/or system of the invention, several selected tasks
could be implemented by hardware, by software or by firmware or by
a combination thereof using an operating system.
[0031] For example, hardware for performing selected tasks
according to embodiments of the invention could be implemented as a
chip or a circuit. As software, selected tasks according to
embodiments of the invention could be implemented as a plurality of
software instructions being executed by a computer using any
suitable operating system. In an exemplary embodiment of the
invention, one or more tasks according to exemplary embodiments of
method and/or system as described herein are performed by a data
processor, such as a computing platform for executing a plurality
of instructions. Optionally, the data processor includes a volatile
memory for storing instructions and/or data and/or a non-volatile
storage, for example, a magnetic hard-disk and/or removable media,
for storing instructions and/or data. Optionally, a network
connection is provided as well. A display and/or a user input
device such as a keyboard or mouse are optionally provided as
well.
BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS
[0032] The patent or application file contains at least one drawing
executed in color. Copies of this patent or patent application
publication with color drawing(s) will be provided by the Office
upon request and payment of the necessary fee.
[0033] Some embodiments of the invention are herein described, by
way of example only, with reference to the accompanying drawings
and images. With specific reference now to the drawings in detail,
it is stressed that the particulars shown are by way of example and
for purposes of illustrative discussion of embodiments of the
invention. In this regard, the description taken with the drawings
makes apparent to those skilled in the art how embodiments of the
invention may be practiced.
[0034] In the drawings:
[0035] FIG. 1 is a flowchart diagram of a method suitable for
processing an image according to various exemplary embodiments of
the present invention;
[0036] FIGS. 2A-B show the effect of an isolated objects removal
procedure, performed according to some embodiments of the present
invention;
[0037] FIG. 3 shows a mapping function which can be applied
according to some embodiments of the present invention for
selective contrast enhancement;
[0038] FIGS. 4A-B are schematic illustrations of a center,
surrounding and remote regions, according to some embodiments of
the present invention;
[0039] FIG. 5 is a schematic illustration of a procedure for
calculating a contrast quantifier, according to some embodiments of
the present invention;
[0040] FIGS. 6A-C show the effect of the contrast quantifiers of
some embodiments of the present invention;
[0041] FIGS. 7A-C show noise analysis of a mammogram, performed
according to some embodiments of the present invention;
[0042] FIG. 8 is a simplified illustration of a system for
processing an image, according to some embodiments of the present
invention;
[0043] FIG. 9 is a simplified illustration of an imaging and
processing system, according to some embodiments of the present
invention;
[0044] FIGS. 10A-B show representative examples of mapping
functions as constructed according to some embodiments of the
present invention for a KODAK mammography system (FIG. 10A) and a
GE mammography system (FIG. 10B);
[0045] FIG. 11 is a flowchart diagram describing a companding
procedure used in experiments performed in accordance with some
embodiments of the present invention;
[0046] FIGS. 12A-F show representative examples of unprocessed
mammograms (FIGS. 12A, 12C and 12E) and the same mammograms
following image processing performed in accordance with some
embodiments of the present invention (FIGS. 12B, 12D and 12F).
[0047] FIGS. 13A-B show Receiver Operating Characteristic (ROC)
analysis applied to data received from an expert radiologist;
and
[0048] FIGS. 14A-B show Receiver Operating Characteristic (ROC)
analysis applied to data received from another expert
radiologist;
DESCRIPTION OF SPECIFIC EMBODIMENTS OF THE INVENTION
[0049] The present invention, in some embodiments thereof, relates
to image processing and, more particularly, but not exclusively, to
a method and system for companding an image.
[0050] Before explaining at least one embodiment of the invention
in detail, it is to be understood that the invention is not
necessarily limited in its application to the details of
construction and the arrangement of the components and/or methods
set forth in the following description and/or illustrated in the
drawings and/or the Examples. The invention is capable of other
embodiments or of being practiced or carried out in various
ways.
[0051] The present embodiments are concerned with method and system
for processing an image to facilitate its display or, in some
embodiments of the present invention, to allow recognizing in the
image features that would not be recognized without the inventive
processing. At least part of the processing can be implemented by a
data processing system, e.g., a computer, configured for receiving
the image and executing the operations described below.
[0052] The method of the present embodiments can be embodied in
many forms. For example, it can be embodied in on a tangible medium
such as a computer for performing the method operations. It can be
embodied on a computer readable medium, comprising computer
readable instructions for carrying out the method operations. In
can also be embodied in electronic device having digital computer
capabilities arranged to run the computer program on the tangible
medium or execute the instruction on a computer readable
medium.
[0053] Computer programs implementing the method of the present
embodiments can commonly be distributed to users on a distribution
medium such as, but not limited to, a floppy disk, a CD-ROM, a
flash memory device and a portable hard drive. From the
distribution medium, the computer programs can be copied to a hard
disk or a similar intermediate storage medium. The computer
programs can be run by loading the computer instructions either
from their distribution medium or their intermediate storage medium
into the execution memory of the computer, configuring the computer
to act in accordance with the method of this invention. All these
operations are well-known to those skilled in the art of computer
systems.
[0054] The image to be analyzed using the teachings of the present
embodiments is generally in the form of imagery data arranged
gridwise in a plurality of picture-elements (e.g., pixels, group of
pixels, etc.).
[0055] The term "pixel" is sometimes abbreviated herein to indicate
a picture-element. However, this is not intended to limit the
meaning of the term "picture-element" which refers to a unit of the
composition of an image.
[0056] References to an "image" herein are, inter alia, references
to values at picture-elements treated collectively as an array.
Thus, the term "image" as used herein also encompasses a
mathematical object which does not necessarily correspond to a
physical object. The original and processed images certainly do
correspond to physical objects which are the scene from which the
imaging data are acquired.
[0057] Each pixel in the image can be associated with a single
digital intensity value, in which case the image is a grayscale
image. Alternatively, each pixel is associated with three or more
digital intensity values sampling the amount of light at three or
more different color channels (e.g., red, green and blue) in which
case the image is a color image. Also contemplated are images in
which each pixel is associated with a mantissa for each color
channels and a common exponent (e.g., the so-called RGBE format).
Such images are known as "high dynamic range" images.
[0058] In some embodiments of the invention, the image is an X-ray
image, also known as a roentgenogram or radiograph. These
embodiments are particularly useful for medical application. For
example, the image can be an X-ray image of a breast of a human
subject, in which case the image is referred to as a mammography
image or a mammogram.
[0059] It is generally accepted that mammography is an effective
and reliable procedure in the early detection of breast cancer.
Mammography is typically performed using X-ray, but may also be
performed using other imaging modalities which are not excluded
from the scope of the present invention. It is recognized by the
present inventors that traditional mammography does not always
provide adequately high-quality images to detect cancer,
particularly in the relatively large population of women having
radiodense breast tissue (e.g., younger women). Unlike other X-ray
images, mammograms require high-quality images because the tissue
density between adipose (fatty), glandular, calcified, or cancerous
tissue is less diverse than, for example, flesh and bone. The
present inventors recognized that subtler contrasts are desirable
to distinguish between these types of tissue.
[0060] Traditional film mammograms have a non-linear response to
x-ray exposure. That is, for example, doubling the x-ray exposure
of film or halving the breast density, does not result in an image
that is twice as bright. As a result, a single traditional film
x-ray exposure often does not show the entire tonal range of a
patient's breast tissue. Often, a radiologist may take exposures at
different energy levels to provide images with different contrasts.
This exposes the patient to several doses of x-rays. Recently,
digital mammography has been proposed as a technology which
replaces the phosphor/film detector with a digital image detector,
with the prospect of overcoming some of the limitations of
film-screen mammography in order to provide higher quality
mammographic images. The advantage of digital mammography is that
it allows separating the image acquisition from the image display.
An additional advantage is that digital detectors provide a much
greater range of contrast than film and the contrast response is
linear over the entire range. This allows digital detectors to more
easily distinguish subtle differences in attenuation of x-rays as
they pass through soft tissue in the breast. Other advantages of
digital mammography include digital image archival and image
transmission to remote locations for viewing purposes.
[0061] Thus, in various exemplary embodiments of the invention the
image is a digital mammogram. Optionally, the processing is
performed such as to allow the radiologist to view the entire
information simultaneously without having to adjust the gain. A
mammogram processed according to this embodiment is advantageous
over traditional mammograms since in traditional digital mammograms
(either mammograms which are acquired digitally, or mammograms that
are acquired analogically and sampled digitally thereafter) the
dynamic range is typically larger than that of the display device.
For example, when a digital mammogram has 2.sup.12=4096 intensity
levels, conventional displays, featuring 2.sup.8-2.sup.10 gray
levels are insufficient for displaying the full dynamic range of
the mammogram, and a substantial portion of the information is lost
once the mammogram is displayed. Thus, in conventional mammograms,
the radiologist oftentimes has to adjust the gain of the display
device such that each time a different portion of the dynamic range
is displayed.
[0062] The present embodiments are also useful for processing other
digital images, particularly images in which the intensity or
grey-level values are relative and not absolute. Thus, besides
mammograms, images suitable for processing according to the
technique of the present embodiments include, but are not limited
to, X-ray images other than CT images, gamma images, PET images,
SPECT images, thermal images and the like.
[0063] Referring now to the drawings, FIG. 1 is a flowchart diagram
of a method suitable for processing an image according to various
exemplary embodiments of the present invention. It is to be
understood that, unless otherwise defined, the operations described
hereinbelow can be executed either contemporaneously or
sequentially in many combinations or orders of execution.
Specifically, the ordering of the flowchart diagrams is not to be
considered as limiting. For example, two or more operations,
appearing in the following description or in the flowchart diagrams
in a particular order, can be executed in a different order (e.g.,
a reverse order) or substantially contemporaneously. Additionally,
several operations described below are optional and may not be
executed.
[0064] The method begins at 10 and continues to 11 at which an
image is received, preferably in a digital form. The method
optionally and preferably continues to 12 at which the image is
processed for removing isolated objects from the image. It is
recognized by the present inventors that isolated elements can
yield a high local and remote contrast responses, which contrasts
may screen other objects from the image. Thus, the removal of
isolated elements moderates the contrasts in the image and
facilitates the companding operation at later stages of the
processing.
[0065] In various exemplary embodiments of the invention the
isolated objects are removed by employing a connected components
labeling procedure. In the labeling procedure each picture-element
is checked to determine if it is connected to adjacent
picture-element, all belonging to the same connected component. If
the picture-element is determined to be connected to an already
identified component, it is labeled with the same component number
given to that component, and if the picture-element is determined
not to be connected to an already identified component, it is given
or labeled with a new component number.
[0066] Once all or most of picture-elements in an image are checked
the largest component or several largest components in the labeled
image is preferably identified and kept, while all other components
are removed from the image. This provides a new image which only
includes the non-removed object or objects. Such image is referred
to herein as segmented image. A connected components labeling
procedure suitable for the present embodiments is found, e.g., in
Gonzalez et al., (2003) Digital image processing using MATLAB,
55-69, 65-76, Prentice-Hall, Inc. Upper Saddle River, N.J.,
USA.
[0067] In various exemplary embodiments of the invention 12 is
preceded by a down-sampling process selected to reduce the number
of identified components in the image. This operation is preferably
employed when it is desired to reduce the processing time. However,
when the computing resources are sufficient, such down-sampling is
not necessary. Also contemplated, are embodiments in which the
input data are normalized, preferably prior to the execution of 12.
The normalization is optional and serves for computation
convenience. Normalization can be performed by linear
transformation which is bases on the principles of a cumulative
distribution function (CDF). Thus, for example, the lowest value of
input intensity can be mapped to a predetermined percentage at the
left hand side of the CDF (e.g., the 5th percentage) and the
highest value of input intensity is mapped to a predetermined
percentage at the right hand side of the CDF (e.g., the 95th
percentage). Formally, such linear transformation can be written
as:
I -> I_high - I_low I_ 0.95 - I_ 0.05 ( I - I_ 0.95 ) + I_high (
EQ . 1 ) ##EQU00001##
where I_high and I_low are, respectively, the highest and lowest
input intensity levels and I.sub.--0.05 and I.sub.--0.95 are,
respectively, the intensity levels of the CDF at the 5th and 95
percentages.
[0068] A representative example of the effect of 12 as applied to a
mammogram is shown in FIGS. 2A and 2B, where FIG. 2A shows the
input mammogram and FIG. 2B shows the largest connected object
after the removal of isolated objects from the mammogram.
[0069] The method optionally and preferably continues to 13 at
which a selective contrast enhancement is applied to enhance the
contrast of some picture-elements of the image. In various
exemplary embodiments of the invention the contrast of
picture-element with low intensity levels is enhanced. In some
embodiments, image data with intensity levels below some intensity
threshold are diminished. Image data with sufficiently high
intensity levels preferably remains unchanged. These embodiments
preferably feature a mapping function having a plateau region below
the intensity threshold.
[0070] A representative example of such mapping function is
illustrated in FIG. 3, showing the dependence of the mapping
function on the intensity levels. As shown, the mapping function
has a plateau for intensity levels below the intensity threshold
and is monotonic for intensity levels above the intensity
threshold. The value of the mapping function at the plateau is such
that once the function is applied to the image data, any
picture-element with intensity level below the intensity threshold
receives a new intensity level which is zero or close to zero
(e.g., below 0.05 or below 0.01) at the plateau region.
[0071] A mathematical expression suitable to be used as a mapping
function is:
I exp ( x , y ) = 1 1 + - a ( I ( x , y ) - c ) I ( x , y ) b , (
EQ . 2 ) ##EQU00002##
where the pair (x,y) is the location of the corresponding
picture-element over the image, I is the input intensity level of
(x, y), I.sub.exp is the new assigned intensity level of (x,y), and
a, b and c are constants. The constant a is associated with the
slope of the curve, the constant c is associated with the intensity
threshold and the constant b is a associated with the curvature of
the curve.
[0072] The intensity threshold can be predetermined or, more
preferably, it can be calculated based on the input image data
before the contrast enhancement 13 is executed.
[0073] In some embodiments of the present invention the intensity
threshold is an average intensity over all picture-elements of the
image having intensity levels within a predetermined intensity
range. These embodiments are particularly useful for mammograms
which are acquired digitally. The predetermined intensity range is
typically expressed as a fraction of the entire range of input
intensities. For example, the predetermined intensity range can
include all intensities between percentage X and percentage Y,
where X and Y are percentage values, e.g., X can be from about 1 to
about 5 and Y can be from about 28 to about 32. When the
aforementioned normalization is employed the predetermined
intensity range can be derived from the CDF of the image.
[0074] A mathematical expression suitable for calculating the
intensity threshold, c, according to some embodiments of the
present invention is:
c = .intg. I s _ 0.02 I s _ 0.3 I s f ( I s ) I s .intg. I s _ 0.02
I s _ 0.3 f ( I s ) I s , ( EQ . 3 A ) ##EQU00003##
where f represents the CDF and I.sub.s represents the input
intensity level. The integration limits I.sub.s.sub.--0.02 and
I.sub.s.sub.--0.3 are intensity levels corresponding to the 2nd and
30th percentages of the CDF. Thus, I.sub.s.sub.--0.02 satisfies
f(I.sub.s.sub.--0.02)=0.02 and I.sub.s.sub.--0.3 satisfies
f(I.sub.s.sub.--0.3)=0.3. Percentages other than 2 and 30 are not
excluded from the scope of the present invention. When Equation 3A
is used for calculating the intensity threshold c, the values of
the parameters a and b in Equation 2 are preferably fixed. Typical
values for these parameters are from about 0.001 to about 0.05 for
a, and from about 0.2 to about 2 for b.
[0075] In some embodiments of the present invention the intensity
threshold is calculated based on the fraction RL of
picture-elements having intensity levels within a predetermined
intensity range. These embodiments are particularly useful for
mammograms which are acquired analogically but are thereafter
converted to digital images. The predetermined intensity range
depends on the type of imaging system. For example, for a KODAK
mammography system, the range can be from about 400 to about 800,
but it is not intended to limit the scope of the present invention
for this range. RL is the number of picture-elements having
intensity levels within the predetermined intensity range as a
fraction of the total number of picture-elements. Once RL is
determined, the intensity threshold is preferably calculated as a
linear function of RL, e.g.,
c=p.sub.0+p.sub.1RL (EQ. 3B)
where p.sub.0 and p.sub.1 are constant parameters. Typical values
for p.sub.0 and p.sub.1 are p.sub.0=1100 and p.sub.1=-100. When
Equation 3B is used for calculating the intensity threshold c, the
values of the parameters a and b in Equation 2 are preferably also
calculated as linear functions of RL. Representative examples of
such linear functions are provided in the Examples section that
follows.
[0076] Also contemplated, are combination of embodiments, for
example, calculation of c using Equation 3A and calculation of a
and/or b using linear functions of RL, or calculation of c using
Equation 3B and fixing the values of a and/or b,
[0077] When the image is an image of a section of a living body the
intensity threshold c is preferably correlated with intensity
levels characterizing fatty tissue present in body section. This
embodiment is particularly useful for mammograms where the fatty
tissue generally has low contrast, and it is desired to selectively
enhance the contrast at picture-elements corresponding to fatty
tissue.
[0078] At 14, a companding procedure is employed so as to reduce
the dynamic range of the image to a second dynamic range which is
displayable on a display device.
[0079] Generally, the companding procedure, shrinks the prominent
gradients over the image (compression) but emphasized small details
therein (expansion). The companding can be according to any
companding procedure known in the art. Representative examples
include the companding procedure described in U.S. Published
Application No. 20040165086, International Patent Publication Nos.
WO2004/075535 and WO2009/081394, the contents of which are hereby
incorporated by reference. Also contemplated is the companding
procedures disclosed in Yuanzhen et al. "Compressing and companding
high dynamic range images with subband architectures,"
International Conference on Computer Graphics and Interactive
Techniques Proceedings of ACM SIGGRAPH 2005 (2005) pages 836-844,
and Duan et al., "Comprehensive Fast Tone Mapping for High Dynamic
Range Image Visualization," Proceedings of Pacific Graphics (2005).
Also contemplated are combinations of the techniques disclosed in
two or more of the above documents. A preferred companding
procedure will now be explained.
[0080] The following procedure is inspired by a known physiological
phenomenon in the human vision system which is called "induction
phenomenon." According to this mechanism, the perceived image is
not just a simple function of the stimulus from specific receptor
cell, but rather a more complicated combination of other stimuli,
originating from other receptors in the field of view.
[0081] Hence, a given ganglion cell in the retina receives input
both from receptors at the center receptive field area and
receptors of nearby, referred to below as surrounding, receptive
field area. The human vision mechanism operates in the retinal
ganglion cell by subtracting surround responses from center
responses. In addition, it is believed that the processing at the
retinal ganglion cell level further includes influences of
responses from receptors being in a "remote" area of the receptive
field that is even farther than the "surround" area from the
"center" area. These influences are typically through the amacrine
cells and horizontal cells of the retina.
[0082] The companding procedure of the present embodiments is
described below for an arbitrarily chosen picture-element,
generally referred to herein as element 40. It is to be understood
that various embodiments of the invention are applied independently
for most or all the picture-elements of the image.
[0083] The companding procedure of the present embodiments
preferably defines several regions in the vicinity of
picture-element 40 (but not necessarily adjacent thereto) are
defined. In various exemplary embodiments of the invention a
surrounding region and a remote region are defined for element 40.
In some embodiments, a center region comprises element 40 and
picture elements immediately adjacent to element 40, is also
defined. Alternatively, the center region may be a single element
region, hence comprising only element 40. This alternative, of
course, coincide with the embodiment in which no center region is
selected.
[0084] As further detailed hereinafter, the regions, preferably
measured in units of picture-elements, are used to simulate the
physiological adaptation mechanism of the human visual system.
[0085] The concept of the center, surrounding and remote regions
may be better understood from the following example, with reference
to FIGS. 4A-B. Thus, if the picture elements are arranged in a
rectangular grid 46, the center region may be a single picture
element (element 40), the surround region may be picture elements
42 surrounding picture elements 40 and the remote region may be
picture elements 44, surrounding picture elements 42.
[0086] In FIG. 4A, the surround region comprises eight
picture-elements immediately surrounding (i.e., adjacent to)
element 40, and the remote region comprises 40 picture-element
forming the two layers surrounding those eight picture-elements.
However, this need not necessarily be the case, since, for some
applications, it may be desired to extend the surround region
farther from those eight elements which immediately surround
element 40. FIG. 4B, for example, illustrates an embodiment in
which the surround region comprises 48 picture-element which form
the first three layers surrounding element 40, and the remote
region comprises 176 picture-element which form the four layers
surrounding those 48 elements. Also contemplated are embodiments in
which the center region comprises more picture-element, e.g., an
arrangement of 2.times.2 or 3.times.3 picture-elements. Other
definitions for the center, surrounding and remote regions are not
excluded from the present invention, both for a rectangular grid or
for any other arrangement according to which the picture elements
of the image are inputted.
[0087] Once the surround and optimally the remote and/or center
regions are defined, the intensities of the picture elements in
each region are preferably used for calculating, for each region,
an overall regional intensity, G.sub.r, where the subscript "r" is
to be understood as a regional subscript. Specifically, for the
center region r should be replaced by the subscript "center", for
the surrounding region r should be replaced by the subscript
"surround" and for the remote region r should be replaced by the
subscript "remote".
[0088] According to a preferred embodiment of the present invention
the overall intensity may be calculated using a regional spatial
profile, f.sub.r. More preferably, the overall intensity is
calculated as an inner product or convolution of the intensity
level I of the picture elements in each region with the respective
regional spatial profile.
[0089] Preferably, the intensity level following the aforementioned
selective contrast enhancement is employed. Mathematically this
inner product or convolution is realized by the following
equation:
G.sub.r=.intg.I*f.sub.rds (EQ. 4)
where the integration is to be understood as convolution, and where
ds is an area integration measure, which is selected in accordance
with the arrangement of the inputted picture elements, e.g., for a
rectangular x-y grid-like arrangement ds equals dx dy. The regional
spatial profiles, f.sub.r, used for calculating the overall
regional intensities are preferably spatial decaying functions,
with may have different forms, depending on the region in which the
profiles are applied. Examples for the specific form of each
regional spatial profile include, but are not limited to, a
Gaussian, an exponent, a Lorenzian, a modified Bessel function and
a power-decaying function. These function are typically
characterized by a slope, denoted K.sub.r (r="center", "surround",
"remote"), which may be different for each region.
[0090] The integration of Equation 4 spans over the respective
region. Typically, but not necessarily, for the overall center
intensity, G.sub.center, the integration spans over a single
picture element (element 40). For example, when a picture-element
is a pixel, G.sub.center can equal the intensity level associated
with this pixel. For the overall surround intensity G.sub.surround,
the integration can extend over a radius of 3 picture-elements, not
including element 40, and for the overall remote intensity,
G.sub.remote, the integration area can extend over an annulus
having an inner radius of 3 picture-elements and an outer radius of
7 picture-elements (see, e.g., FIG. 4B). At the boundaries of the
images, all the integrations preferably facilitate periodic
boundary conditions.
[0091] In various exemplary embodiments of the invention one or
more adaptation terms are calculated using intensity levels of each
picture element of the surrounding and remote regions. These
adaptation terms are denoted herein as .sigma..sub.r, where r is a
regional index, r, as described above. In various exemplary
embodiments of the invention two adaptation terms are defined,
.sigma..sub.center and .sigma..sub.surround.
[0092] Each of the adaptation terms .sigma..sub.center and
.sigma..sub.surround preferably comprises a local part and a remote
part, to account for the influence of the remote region:
.sigma..sub.r=.sigma..sub.r,local+.sigma..sub.r,remote, (EQ. 5)
where "r"="center" or "surround". .sigma..sub.center, local is
preferably a function of the overall center intensity G.sub.center,
and .sigma..sub.center, remote is preferably a function of overall
remote intensity G.sub.remote. Similarly, .sigma..sub.surround,
local is preferably a function of the overall surround intensity
G.sub.surround, and .sigma..sub.surround, remote is preferably a
function of overall remote intensity G.sub.remote. The local parts
and the remote parts of the center and surround adaptation terms
determine the relative weight of the remote region on each
adaptation term.
[0093] The local parts of .sigma..sub.center and
.sigma..sub.surround are preferably calculated as linear functions
of the regional overall intensities (with either constant or
variable coefficients).
.sigma..sub.r,local=.alpha..sub.rG.sub.r+.beta..sub.r, (EQ. 6)
The value of the coefficients .alpha..sub.r and .beta..sub.r
("r"="center" or "surround") may be selected in accordance with the
dynamic range of the image. The coefficient .beta..sub.r generally
controls the amount of adaptation which is performed by the
adaptation terms. More specifically, as can be understood by
inspecting Equation 6, .beta..sub.r functions as an adaptation
threshold, where for large values of .beta..sub.r, the adaptation
level is low and for low values of .beta..sub.r, the adaptation
level is high. The amount of the required adaptation can also
depend on the intensity at a particular picture element or on the
intensity of a group of picture elements, e.g., the picture
elements populating the remote region or more. Thus, in these
embodiments .beta..sub.r is a function rather than a constant
coefficient. According to some of these embodiments, .beta..sub.r
is an increasing function of the intensity, so that for low
intensity regions the amount of adaptation is high and for high
intensity regions the amount of adaptation is low.
[0094] Although the coefficients of Equation 6 are mathematical
coefficients, they are preferably based on electro-physiological
findings. In accordance with the above physiological "gain
control", each of the center and surround adaptation term
independently characterizes a dynamical intensity curve. The
coefficients .alpha..sub.r are thus determine the degree of
curve-shifting, for example, higher values of .alpha..sub.r lead to
higher shifting amount of the response curve. The combination
between .alpha..sub.r and .beta..sub.r determine the onset of the
gain-control mechanism, hence .alpha..sub.r and .beta..sub.r serve
as gain thresholds.
[0095] The remote parts of .sigma..sub.center and
.sigma..sub.surround can be calculated using one or more adaptation
weights denoted MRM.sub.r ("r"="center" or "surround"):
.sigma..sub.r,remote=MRM.sub.rG.sub.remote. (EQ. 7)
[0096] Generally, MRM.sub.r can be a constant or it can vary with
one or more of the regional intensities. In some embodiments of the
present invention the same adaptation weight MRM.sub.r is used for
the calculation of the remote parts of both .sigma..sub.center and
.sigma..sub.surround, and in some embodiments of the present
invention the adaptation weight used for the calculation of the
remote part of .sigma..sub.center differ from that used for the
calculation of the remote part of .sigma..sub.surround. For
example, in some embodiments the adaptation weight used for the
calculation of the remote part of .sigma..sub.center is a function
of one or more of the regional intensities and the adaptation
weight used for the calculation of the remote part of
.sigma..sub.surround is a constant.
[0097] As can be seen form Equation 7, the remote parts of
.sigma..sub.center and .sigma..sub.surround are defined as a
multiplication between the adaptation weight and the overall remote
intensity. The adaptation weights modulate the adaptation in
accordance with the intensity levels of picture elements in the
remote region. A particular feature of the adaptation weights is
that these functions serve also for preserving and/or improving the
contrast of the image, in a manner that will be now described.
[0098] At the vicinity of element 40, a local-contrast can be
defined as the difference between the intensity level of element 40
and the picture elements of the surrounding region, where a large
difference is interpreted as a higher local-contrast and a low
difference is interpreted as a lower local-contrast. This
difference may be calculated by more than one way, for example, by
subtraction, division, subtraction of logarithms, discrete
differentiation, discrete logarithmic differentiation or any other
suitable mathematical operation between intensity levels. The
adaptation weights can be functions which are selected in
accordance with the local-contrast calculations. Specifically, in
regions in which the local-contrast is high these functions have
small numerical values and in regions in which the local-contrast
is low these functions have higher numerical values.
[0099] The adaptation weight of the present embodiments describes,
at least in part, contrast induction between the surrounding region
and the remote region of element 40. Optionally and preferably the
adaptation weight features a plurality of different resolutions,
wherein a contrast induction is calculated separately for each
resolution. For example, for each resolution, contrast inductions
are calculated using subtractions between intensity levels
integrated over the surrounding region and intensity levels
integrated over the central region.
[0100] Mathematically, this procedure can be expressed as follows.
Let L.sup.k.sub.r be an intensity integrated over the region r
(r=center or surround) at resolution k. A representative expression
for L.sup.k.sub.r can be:
L.sup.k.sub.r=.intg.I*f.sub.r.sup.kds, (EQ. 8)
where the integration of spans over the respective (center or
surround) region, f.sub.r.sup.k are spatial profiles at resolution
k, which may have the same properties as the spatial profiles
f.sub.r of Equation 4, and I is the intensity level of element 40,
preferably following the selective contrast enhancement procedure
described above. Thus, in various exemplary embodiments of the
invention I=I.sub.exp. The resolution of the spatial profiles is
preferably expressed as a local support of the respective profile.
For example, when the profile is a Gaussian profile, f.sub.r.sup.k
can be written as (r=center or surround):
f r k ( x , y ) = 1 .pi. .rho. r k exp ( - x 2 + y 2 ( .rho. r k )
2 ) ( EQ . 9 ) ##EQU00004##
where .rho..sub.r.sup.k is the local support of region r for
resolution k. Typically, the surround profile for a given
resolution has a three times larger support area and a radius which
is two times larger than of that of the center profile. Thus, the
present embodiments calculate two series of integrated intensities
L.sup.k.sub.center (k=1, 2, . . . ) and L.sup.k.sub.surround (k=1,
2, . . . ). For each series, each element features a profile of
different support, and therefore each such series is based on a
series of profiles which is based on a series of radii.
Representative radii examples of the series of radii
.rho..sup.k.sub.center and .rho..sup.k.sub.surround are provided in
the Examples section that follows.
[0101] Once the integrated intensities L.sup.k.sub.r are
calculated, they are preferably subtracted to provide a contrast
induction L.sup.k.sub.tot for each resolution k:
L.sup.k.sub.tot=L.sup.k.sub.center-L.sup.k.sub.surround (EQ.
10)
[0102] The contrast inductions are then used for calculating a
local contrast quantifier C.sub.local and a remote contrast
quantifier C.sub.remote. C.sub.local can be calculated as a sum
over all resolution of an inner product or convolution of the
absolute value of L.sup.k.sub.tot with a local spatial profile
W.sub.local:
C local ( i , j ) = local area ( k = 1 num of scales L total k ) (
x , y ) W local ( i - x , j - y ) . ( EQ . 11 ) ##EQU00005##
[0103] C.sub.remote can be calculated as a sum over all resolution
of an inner product or convolution of C.sub.local with a remote
spatial profile W.sub.remote:
C remote ( i , j ) = remote area C local ( x , y ) W remote ( i - x
, j - y ) ( EQ . 12 ) ##EQU00006##
Note that in Equations 11 and 12, the convolutions are expressed in
a discrete form. The skilled person provided with the details
described herein would know how to adjust the convolution to a more
general case by replacing the summation over the local area with an
integration as further detailed hereinabove. The local and remote
spatial profile can have properties are preferably which are
similar to the other profiles described above, but their support
area is preferably larger. In some embodiments of the present
invention W.sub.local has a Gaussian shape and W.sub.remote has an
annular Gaussians shape:
W local ( x , y ) = 1 .pi. .rho. local exp ( - x 2 + y 2 .rho.
local 2 ) W remote ( x , y ) = 1 .pi. .rho. remote exp ( - x 2 + y
2 .rho. remote 2 ) . ( EQ . 13 ) ##EQU00007##
[0104] The radius .rho..sub.local is preferably larger, e.g., two
times larger than .rho..sup.k.sub.center for the coarsest
resolution, and the radius .rho..sub.remote is preferably larger,
e.g., two times larger than .rho..sup.k.sub.surround for the
coarsest resolution.
[0105] Once the local and remote contrast inductions are
calculated, the center and surround adaptation weights can be
calculated as follows:
MRM.sub.center=C.sub.local,max-C.sub.local
MRM.sub.remote=C.sub.remote,max-C.sub.remote (EQ. 14)
where C.sub.local, max and C.sub.remote, max are the maximal values
of C.sub.local and C.sub.remote over a region of predetermined
size. A typical radius of such region is about two- or three times
the radius of the remote region.
[0106] A visual representation of the contrast quantifier
calculation is schematically illustrated in FIG. 5. The outer large
circle represents the "remote" area that is defined by the Gaussian
profile W.sub.remote, the inner circle represents the "local" area
that is defined by the Gaussian profile W.sub.local, and the
distributed small circles of different sizes illustrate the
"center-surround" profiles of different resolutions.
[0107] The effect of contrast quantifiers for a mammogram is shown
in FIGS. 6A-C, where FIG. 6A shows the input mammogram, FIG. 6B
shows the local contrast quantifier C.sub.local(x,y), and FIG. 6C
shows the remote contrast quantifier C.sub.remote(X,Y).
[0108] Once the adaptation terms .sigma..sub.center and
.sigma..sub.surround are calculated, the companding procedure
preferably assigns a new intensity level to element 40 based the
adaptation terms. The new intensity level is in a second dynamic
range which is smaller than dynamic range of the original image.
The second dynamic range is preferably displayable on a display
device, as further detailed hereinabove.
[0109] The assignment of a new intensity level I.sub.res is
preferably according to the following equation:
I res = d cen G cen G cen + .sigma. cen - d srnd G srnd G srnd +
.sigma. srnd ( EQ . 15 ) ##EQU00008##
where, d.sub.cen and d.sub.srnd are constants.
[0110] In various exemplary embodiments of the invention the method
proceeds to 15 at which a noise reduction procedure is employed.
Generally, the noise reduction serves for reducing excessive
grainess which is that is inconvenient for an observing eye. The
present inventors found that the excessive grainess is originated
from white noise with some localized distribution (e.g., Gaussian
distribution) that existed earlier in the original image but was
not prominent due to high dynamic range of the image. FIGS. 7A-C
illustrate the noise in a mammogram following the application of
operations 11-14, where FIG. 7A is the processed mammogram, FIG. 7B
is an enlarge view of a generally empty area (marked at the lower
right corner of FIG. 7A) and FIG. 7C is a histogram of intensities
over the area of FIG. 7B. As shown, the noise has a localized
distribution over the image. It was found by the present inventor
that an anisotropic diffusion procedure can be employed for
reducing the noise in the image.
[0111] This is particularly useful for a mammogram since in a
mammogram traditional noise reduction techniques, such as moving
average, frequency domain LPF and Wiener filter may result it loss
of valuable information for correct diagnosis. A suitable
anisotropic diffusion procedure for the present embodiments is
found in Mayo et al., 2006, EMBS '06, 28th Annual International
Conference of the IEEE, and Weickert, J., 2001, Pattern
Recognition, 34 (9), 1813-1824.
[0112] An anisotropic diffusion procedure suitable for the present
embodiments will now be described.
[0113] The model considered herein for noisy image is defined
as:
f={tilde over (f)}+n
where f is the observed image, and it is a sum of an ideal
noise-free image {tilde over (f)} and a noise signal n. It is
assumed that noise is a white noise with Gaussian distribution
centered at zero.
[0114] A variational method for image restoration from noisy images
is obtained as a minimizer of energy functional:
s [ u ] = .intg. .OMEGA. [ .PHI. ( .gradient. u ) + 1 2 .mu. ( u -
f ) 2 ] x y ##EQU00009##
where u is the filtered image and .mu. is a Lagrange multiplier.
The first summand rewards smoothness, and the second summand
encourages similarity between the restored image and the original
one. The second summand is derived from the following
constraint:
.intg. .OMEGA. ( u - f ) 2 x y = .sigma. 2 .intg. .OMEGA. x y
##EQU00010##
where .sigma. is the standard deviation of the noise. The regulizer
.PHI. is defined as:
.PHI.(|.gradient.u|)= {square root over
(.beta..sup.2+|.gradient.u|.sup.2)}+.epsilon.|.gradient.u|.sup.2
where .beta. and .epsilon. are constant parameters. The first
summand of .PHI. is responsible for anisotropic diffusion, and the
second summand is for linear diffusion.
[0115] From variational calculus, it follows that the minimizer of
s[u] satisfies the Euler-Lagrange equation:
div ( .PHI. ' ( .gradient. u ) .gradient. u .gradient. u ) + .mu. (
f - u ) = 0 ##EQU00011##
[0116] This can be regarded as the steady-state (t.fwdarw..infin.)
of the diffusion process given by gradient descend:
u t = - s [ u ] u = div ( .PHI. ' ( .gradient. u ) .gradient. u
.gradient. u ) + .mu. ( f - u ) ##EQU00012##
where
.PHI. ' ( .gradient. u ) .gradient. u ##EQU00013##
represents a diffusion coefficient. .PHI. derivative is:
.PHI. ' ( .gradient. u ) = .gradient. u .beta. 2 + .gradient. u 2 +
.gradient. u ##EQU00014##
[0117] Thus, the total anisotropic diffusive process with Newman
boundary condition is:
{ u t = div ( .gradient. u .beta. 2 + .gradient. u 2 ) + .gradient.
2 u + .mu. ( f - u ) u n | .differential. .OMEGA. = 0
##EQU00015##
[0118] Full development of diffusive process equations and its
discretization is provided in the Examples section that
follows.
[0119] The stopping criterion for the diffusive process can be as
follows. Assuming that the unknown noise n is uncorrelated with the
unknown signal {tilde over (f)}, the noise reduction procedure
minimizes the covariance of the "noise" f-u(t) with the filtered
image u(t), or employ its normalized form, the correlation
coefficient:
corr(f-u(t),u(t))
[0120] The stopping time T can be selected such that the expression
is as small as possible. Thus, In some embodiments of the present
invention the stopping criterion is defined as:
T = arg min t corr ( f - u ( t ) , u ( t ) ) ##EQU00016##
[0121] At 16 the method outputs the processed image to a computer
readable medium or a display device.
[0122] The method ends at 17.
[0123] Reference is now made to FIG. 8, which is a simplified
illustration of a system 70 for processing an image, according to
various exemplary embodiments of the present invention. In some
embodiments of the present invention system 70 comprises an input
unit 72 for inputting a plurality of picture elements which form
the image. As further detailed hereinabove, the plurality of
picture elements contains intensity values, and includes the
arbitrarily chosen element 40, to which the following description
refers. System 70 comprises a preprocessing unit 74 which
preprocesses the input image to remove isolated objects from the
image as further detailed hereinabove. In various exemplary
embodiments of the invention unit 74 is configured for executing
one or more of the operations of method 10 described above.
[0124] System 70 further comprises a companding unit 78 which
applies a companding procedure to the preprocessed image, so as to
reduce the dynamic range of the image. In various exemplary
embodiments of the invention unit 78 is configured for executing
one or more of the operations of the companding procedure described
above. In some embodiments of the present invention the dynamic
range is reduced to a second dynamic range which is displayable on
a display device 80. In some embodiments of the present invention
system 70 comprises display device 80.
[0125] Reference is now made to FIG. 9 which is a simplified
illustration of an imaging and processing system 90, according to
various exemplary embodiments of the present invention. System 90
preferably comprises an imaging system 92 which acquire an image of
a living subject. Preferably, system 92 acquires a mammogram, but
other types of images are not excluded from the scope of the
present invention. Preferably system 92 acquires digital images in
which the intensity or grey-level values are relative and not
absolute. Representative examples for system 92 including, without
limitation, a mammography system, an X-ray system other than a CT
scanner, a gamma camera system, a PET system, a SPECT system, a
thermal imaging system and the like.
[0126] System 90 further comprises an image processing system 94
which processes the image and provides a processed image having a
second dynamic range which is smaller than the first dynamic range.
In various exemplary embodiments of the invention image processing
system 94 comprises system 70. In some embodiments image processing
apparatus 94 is system 70.
[0127] As used herein the term "about" refers to .+-.10%.
[0128] The word "exemplary" is used herein to mean "serving as an
example, instance or illustration." Any embodiment described as
"exemplary" is not necessarily to be construed as preferred or
advantageous over other embodiments and/or to exclude the
incorporation of features from other embodiments.
[0129] The word "optionally" is used herein to mean "is provided in
some embodiments and not provided in other embodiments." Any
particular embodiment of the invention may include a plurality of
"optional" features unless such features conflict.
[0130] The terms "comprises", "comprising", "includes",
"including", "having" and their conjugates mean "including but not
limited to".
[0131] The term "consisting of means "including and limited
to".
[0132] The term "consisting essentially of" means that the
composition, method or structure may include additional
ingredients, steps and/or parts, but only if the additional
ingredients, steps and/or parts do not materially alter the basic
and novel characteristics of the claimed composition, method or
structure.
[0133] As used herein, the singular form "a", "an" and "the"
include plural references unless the context clearly dictates
otherwise. For example, the term "a compound" or "at least one
compound" may include a plurality of compounds, including mixtures
thereof.
[0134] Throughout this application, various embodiments of this
invention may be presented in a range format. It should be
understood that the description in range format is merely for
convenience and brevity and should not be construed as an
inflexible limitation on the scope of the invention. Accordingly,
the description of a range should be considered to have
specifically disclosed all the possible subranges as well as
individual numerical values within that range. For example,
description of a range such as from 1 to 6 should be considered to
have specifically disclosed subranges such as from 1 to 3, from 1
to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6 etc., as
well as individual numbers within that range, for example, 1, 2, 3,
4, 5, and 6. This applies regardless of the breadth of the
range.
[0135] Whenever a numerical range is indicated herein, it is meant
to include any cited numeral (fractional or integral) within the
indicated range. The phrases "ranging/ranges between" a first
indicate number and a second indicate number and "ranging/ranges
from" a first indicate number "to" a second indicate number are
used herein interchangeably and are meant to include the first and
second indicated numbers and all the fractional and integral
numerals therebetween.
[0136] It is appreciated that certain features of the invention,
which are, for clarity, described in the context of separate
embodiments, may also be provided in combination in a single
embodiment. Conversely, various features of the invention, which
are, for brevity, described in the context of a single embodiment,
may also be provided separately or in any suitable subcombination
or as suitable in any other described embodiment of the invention.
Certain features described in the context of various embodiments
are not to be considered essential features of those embodiments,
unless the embodiment is inoperative without those elements.
[0137] Various embodiments and aspects of the present invention as
delineated hereinabove and as claimed in the claims section below
find experimental support in the following examples.
EXAMPLES
[0138] Reference is now made to the following examples, which
together with the above descriptions illustrate some embodiments of
the invention in a non limiting fashion.
Example 1
Processing of Mammograms
[0139] Embodiments of the present invention have been utilized for
processing mammograms.
[0140] Images were acquired by KODAK and GE workstations. The KODAK
CR mammography station acquires the images analogically and
converts them to digital, while the GE DR station acquires images
digitally.
[0141] The images were supplied by department of diagnostic imaging
from the Sheba medical center, Israel. The common format for
digital mammograms is the DICOM format. This format contains an
uncompressed image data and a file header which contains
information about the patient, the examination and different image
properties. The spatial resolution of KODAK images used in the
present example is 49 .mu.m/pixel and the matrix size is
4784.times.3584; the matching properties of GE images are 94
.mu.m/pixel and 2294.times.1914. The contrast resolution of both
modalities is 12 bits/pixel (2.sup.12 gray levels). The images were
stored in BMP format.
[0142] The mapping function parameters employed were as follows.
For the GE images a and b were fixed (a=0.011, b=0.8) and c was
calculated adaptively according to Equation 3A above. For the KODAK
images the a, b and c parameters were calculated adaptively based
on the original image histogram. Specifically, RL was the fraction
of pixels with intensity range of 400-800, c was calculated using
Equation 3B with a predetermined intensity range of 400-800,
p.sub.0=1100 and p.sub.1=-100, a was calculated using the linear
relation a=0.04-0.001(RL-1)/6 and b was calculated using the linear
relation b=0.6-(0.8-0.6)(RL-10)/9.
[0143] Representative examples of the mapping function for the
KODAK and GE systems are shown in FIGS. 10A and 10B, respectively.
The normalization parameters for the images were: I_low=750 and
I_high=1250.
[0144] The companding procedure is illustrated graphically in FIG.
11.
[0145] Each pixel of the input image was simulated as the center
region of the opponent receptive field. Accordingly, the center
region had a radius of 1 pixel. The surround region was defined as
having an outer radius of 9 pixels without including the central
pixel. The remote region had an annular shape with inner and outer
radius of 9 and 21 pixels, respectively. The inner radius of the
remote region was equal to the external diameter of the surround
region and therefore did not overlap the center or the surround
regions.
The spatial profiles were Gaussian decaying profiles with
appropriate variances. The total weights of these functions were
normalized to 1. These properties are summarized in Table 1.
TABLE-US-00001 TABLE 1 Region inner radius external radius variance
center -- 1 -- surround 1 9 10 remote 9 21 22
[0146] The parameters for the adaptation terms are summarized in
Table 2, below. These values are the optimal values we could select
over a large number of simulations on a vast variety of mammograms.
The same values were used for all processed mammograms.
TABLE-US-00002 TABLE 2 Parameter Value Parameter Value
.alpha..sub.cen 0.5 .alpha..sub.srnd 0.3 .beta..sub.cen 1.5
.beta..sub.srnd 2 c.sub.cen 1.5 c.sub.srnd 2 d.sub.cen 2.5
d.sub.srnd -0.5
[0147] Local Multi-scale contrast was calculated by 6 resolutions.
The profile's size and its radius increment between two subsequent
resolutions were of 2 pixels. For each resolution the center radius
was of the same size as the profile's size. The radius for the
matching surround profile was twice larger. The different sizes and
the compatible radii at the different resolutions, are summarized
in Table 3.
TABLE-US-00003 TABLE 3 resolution f.sub.cen.sup.k size
.rho..sub.cen.sup.k size f.sub.srnd.sup.k size .rho..sub.srnd.sup.k
size index [pixels] [pixels] [pixels] [pixels] 1 1 1 3 2 2 3 3 9 6
3 5 5 15 10 4 7 7 21 14 5 9 9 27 18 6 11 11 33 22
[0148] The calculations of C.sub.local were performed with a
.rho..sub.local radius (variance of Gaussian profile) of 4 pixels.
The calculations of C.sub.remote quantifier were performed with a
.rho..sub.remote radius of 12 pixels. The sizes of the profiles
W.sub.local and W.sub.remote were 20 and 60 pixels, respectively,
where W.sub.remote has an annular shape with internal radius of 20
pixels.
[0149] For the anisotropic diffusion filter, the following
parameters were used .beta.=1, .epsilon.=0.04, and .DELTA.t=0.1.
The parameter .mu. was calculated iteratively as explained in
Example 3, below.
Example 2
Clinical Test
[0150] The processing method described in Example 1 was applied to
mammograms which were thereafter analyzed by expert
radiologists.
Methods
[0151] The study included 40 female subjects. For each subject 2 or
4 mammograms were acquired (2 of right breast, 2 of left breast or
2 of each breast). The two mammograms of the same breast were the
craniocaudal and the mediolateral-oblique views, both used for the
diagnosis of a single breast.
[0152] 70 different mammograms were chosen for analysis, 20
cancerous mammograms and 50 benign mammograms. The gold standard
for cancerous mammograms was biopsy. The gold standard for benign
mammograms was a two-year follow up of benign results. The
different cases were selected to allow variety of lesions. A high
resolution monitor of 5MP was used to perform the test.
[0153] Two expert radiologists were independently presented with
the mammograms, blinded to patient's identity. Initially, each
expert radiologist was provided with mammograms processed in
accordance with some embodiments of the present invention as
described in Example 1 above. Each expert radiologist analyzed the
processed mammograms and provided an expert diagnosis which
included specifying all lesions found and ranking the respective
mammograms according to the Breast Imaging Reporting and Data
System (BIRADS) scale. The BIRADS scale is presented in Table 4,
below. Expert No. 1 analyzed mammograms Nos. 1-70 and Expert No. 2
analyzed mammograms Nos. 23-70.
[0154] At least two weeks later, the expert radiologists were
provided with the unprocessed versions of the mammograms (namely
mammograms as provided by the mammography systems) for a second
analysis. It is assumed that the second analysis was unbiased by
the results of the first analysis. Following the second analysis
the expert radiologists provided an expert diagnosis which included
specifying all lesions found and ranking the respective mammograms
according to the BIRADS scale. Same as before expert No. 1 analyzed
mammograms Nos. 1-70 and Expert No. 2 analyzed mammograms Nos.
23-70.
The input from the experts was subjected to a Receiver Operating
Characteristic (ROC) analysis [Park et al., 2004, Korean J Radiol,
5(1):11-8]. ROC analysis is suitable for the present study since
the BIRADS scale includes five rating categories. An empirical ROC
curve was obtained by calculating sensitivity and specificity for
four different thresholds .gtoreq.2, .gtoreq.3, .gtoreq.4,
.gtoreq.5. An estimation of smooth (fitted data) ROC curve was made
using maximum likelihood estimation (MLE). To this end, a binormal
distribution was assumed for characterizing the shape of the ROC
curve using two parameters, a=(.mu..sub.1-.mu..sub.0)/.sigma..sub.1
and b=.sigma..sub.0/.sigma..sub.1, where .mu..sub.0, .sigma..sub.0
and .sigma..sub.1.mu..sub.1, .sigma..sub.1 are mean and standard
deviation of the test results without and with lesions,
respectively.
TABLE-US-00004 TABLE 4 BIRADS Rank Finding 0 need additional
imaging evaluation 1 negative 2 benign 3 probably benign - short
interval follow-up suggested 4 suspicious abnormality - biopsy
should be considered 5 highly suggestive of malignancy -
appropriate action should be taken
Results
[0155] Representative examples of some of the mammograms are shown
in FIGS. 12E-F, where FIGS. 12A, 12C and 12E show the unprocessed
mammograms (as provided by the conventional mammography system),
and FIGS. 12B, 12D and 12F show the mammograms after processing in
accordance with some embodiments of the present invention.
[0156] Both experts reported that the processed mammograms
successfully show all biological features that are required for the
diagnosis. In particular, Cooper's ligaments were identified more
easily in the processed mammograms.
[0157] The BIRADS ranks as provided by the experts are summarized
in Table 5. Expert No. 1 tested 70 cases, 20 of which were
cancerous and the rest were benign, and expert No. 2 tested 48
cases, 16 of which were cancerous and the rest were benign.
TABLE-US-00005 TABLE 5 Processed Unprocessed Mammogram Mammogram
Expert Expert Radiologist Radiologist Gold No. 1 No. 2 1 2 Case #
Standard BIRADS BIRADS BIRADS BIRADS 1 benign 3 2 2 benign 1 2 3
benign 2 2 4 benign 2 2 5 benign 3 3 6 benign 1 2 7 benign 3 2 8
benign 3 3 9 benign 3 3 10 benign 3 2 11 benign 1 1 12 benign 1 1
13 benign 3 2 14 benign 3 2 15 benign 1 2 16 benign 3 1 17 benign 3
2 18 benign 2 2 19 cancerous 1 1 20 cancerous 4 4 21 benign 1 1 22
cancerous 2 2 23 benign 2 5 2 4 24 cancerous 5 4 5 4 25 benign 4 2
4 2 26 cancerous 4 4 5 3 27 benign 1 2 1 2 28 cancerous 3 5 2 5 29
benign 1 2 1 2 30 cancerous 5 4 5 3 31 benign 2 5 2 5 32 cancerous
4 2 4 1 33 cancerous 5 2 5 1 34 benign 1 2 2 4 35 benign 1 3 1 4 36
benign 2 1 1 2 37 benign 2 4 4 2 38 benign 1 1 1 2 39 benign 2 1 4
2 40 benign 1 1 1 2 41 benign 1 2 3 2 42 benign 1 3 3 3 43 benign 2
1 2 2 44 benign 2 2 2 2 45 benign 1 2 1 2 46 benign 1 2 1 1 47
benign 1 2 3 1 48 benign 1 4 3 4 49 benign 1 2 1 2 50 benign 2 2 4
2 51 benign 1 3 1 4 52 benign 2 2 1 2 53 benign 2 2 2 2 54 benign 2
3 2 2 55 benign 2 5 2 5 56 benign 2 2 1 3 57 cancerous 5 5 5 5 58
benign 2 5 2 5 59 cancerous 5 4 5 4 60 cancerous 5 3 5 3 61
cancerous 4 5 5 4 62 benign 2 2 2 2 63 cancerous 3 4 4 4 64 benign
2 2 2 2 65 cancerous 5 4 4 4 66 benign 1 3 3 3 67 cancerous 4 4 5 5
68 cancerous 3 5 2 5 69 cancerous 4 4 3 3 70 cancerous 5 3 5 2
[0158] The ROC analyses for Expert Nos. 1 and 2 are shown in FIGS.
13A-B (Expert No. 1) and 14A-B (Expert No. 2). Shown in FIGS. 13A
and 14A are the true positives fraction (TPF) as a function of the
false positives fraction (FPF) obtained for the empirical data.
Shown in FIGS. 13B and 14B are the TPF as a function of the FPF
obtained for the fitted data.
[0159] Tables 6 and 7 summarize the ROC areas as calculated for the
empirical and fitted data for Expert No. 1 (Table 6) and Expert No.
2 (Table 7).
TABLE-US-00006 TABLE 6 Empirical ROC Area Fitted ROC area
unprocessed 0.853 0.867 mammograms processed 0.904 0.913
mammograms
TABLE-US-00007 TABLE 7 Empirical ROC Area Fitted ROC area
unprocessed 0.925 0.937 mammograms processed 0.95 0.958
mammograms
[0160] Tables 8 and 9 summarize the empirical values of the TPF and
FPF as obtained for Expert No. 1 (Table 8) and Expert No. 2 (Table
9), for a threshold of .gtoreq.4 (meaning that BIRADS score which
is 4 or above indicates the presence of cancer).
TABLE-US-00008 TABLE 8 Empirical TPF Empirical FPF unprocessed 0.75
0.22 mammograms processed 0.85 0.18 mammograms
TABLE-US-00009 TABLE 9 Empirical TPF Empirical FPF unprocessed 1
0.25 mammograms processed 1 0.2812 mammograms
[0161] The results presented in this example demonstrate that the
processing technique of the present embodiments increase the
sensitivity without decreasing the specificity. The increment in
ROC curve area demonstrates the effectiveness of the inventive
processing technique.
Example 3
Nonlinear Diffusive Filter
[0162] The model considered herein for noisy image is:
f={acute over (f)}+n
where f is the observed image, and it is a sum of an ideal
noise-free image {tilde over (f)} and a noise signal n. Variational
method is employed for image restoration from noisy images obtain a
filtered version of some degraded image f as the minimizer of
energy functional:
s [ u ] = .intg. .OMEGA. [ .PHI. ( .gradient. u ) + 1 2 .mu. ( u -
f ) 2 ] x y ##EQU00017##
where u is the filtered image and .mu. the Lagrange multiplier. The
first summand rewards smoothness, while the second summand
encourages similarity between the restored image and the original
one. The second summand is derived from the constraint:
.intg. .OMEGA. ( u - f ) 2 x y = .sigma. 2 .intg. .OMEGA. x y
##EQU00018##
where .sigma. is the noise standard deviation. The regulizer .PHI.
is a potential which is defined as:
.PHI.(|.gradient.u|)= {square root over
(.beta..sup.2+|.gradient.u|.sup.2)}+.epsilon.|.gradient.u|.sup.2
where .beta. and .epsilon. are constant parameters. The first
summand of .PHI. is responsible for anisotropic diffusion, and the
second summand is for linear diffusion. From variational
calculus:
ds = lim .fwdarw. 0 s [ u + .eta. ( x , y ) ] - s [ u ] = = .intg.
.OMEGA. [ ( .PHI. ( .gradient. u ) u + 1 2 .mu. ( u - f ) 2 u )
.eta. ( x , y ) + .PHI. ( .gradient. u ) u x .eta. x ( x , y ) +
.PHI. ( .gradient. u ) u y .eta. y ( x , y ) ] x y = = .intg.
.OMEGA. [ .PHI. ( .gradient. u ) u .eta. ( x , y ) + .PHI. (
.gradient. u ) u x .eta. x ( x , y ) + .PHI. ( .gradient. u ) u y
.eta. y ( x , y ) ] x y + .intg. .OMEGA. .mu. ( u - f ) .eta. x y
##EQU00019## .eta. ( x , y ) is a function by which u changes , and
it fulfils .eta. ( x , y ) | .differential. .OMEGA. = 0.
##EQU00019.2##
[0163] The first summand is equal to zero. Therefore:
ds = .intg. .OMEGA. [ .PHI. ' .gradient. u u x .eta. x ( x , y ) +
.PHI. ' .gradient. u u y .eta. y ( x , y ) ] x y + .intg. .OMEGA.
.mu. ( u - f ) .eta. x y = = .intg. .OMEGA. [ .PHI. ' u x
.gradient. u .eta. x ( x , y ) + .PHI. ' u y .gradient. u .eta. y (
x , y ) ] x y + .intg. .OMEGA. .mu. ( u - f ) .eta. x y = = .intg.
.OMEGA. .PHI. ' .gradient. u .gradient. .fwdarw. u .gradient.
.fwdarw. .eta. x y + .intg. .OMEGA. .mu. ( u - f ) .eta. x y
##EQU00020##
[0164] Using Green's first identity for the first summand one
obtains:
ds = .intg. .differential. .OMEGA. .PHI. ' .gradient. u .gradient.
.fwdarw. u ( .eta. n ^ ) - .intg. .OMEGA. div ( .PHI. ' .gradient.
u .gradient. .fwdarw. u ) .eta. x y + .intg. .OMEGA. .mu. ( u - f )
.eta. x y = = .intg. .differential. .OMEGA. .PHI. ' .gradient. u
.eta. ( .gradient. .fwdarw. u n ^ ) - .intg. .OMEGA. div ( .PHI. '
.gradient. u .gradient. .fwdarw. u ) .eta. x y + .intg. .OMEGA.
.mu. ( u - f ) .eta. x y = = .intg. .differential. .OMEGA. .PHI. '
.gradient. u .eta. ( u n .fwdarw. ) - .intg. .OMEGA. div ( .PHI. '
.gradient. u .gradient. .fwdarw. u ) .eta. x y + .intg. .OMEGA.
.mu. ( u - f ) .eta. x y ##EQU00021##
[0165] Using Derichle boundary condition
.eta.(x,y)|.sub..differential..OMEGA.=0. Using Newman boundary
condition
u n .fwdarw. | .differential. .OMEGA. = 0. ##EQU00022##
In any case, the boundary integral equals zero. Therefore:
ds = - .intg. .OMEGA. div ( .PHI. ' .gradient. u .gradient.
.fwdarw. u ) .eta. x y - .intg. .OMEGA. .mu. ( f - u ) .eta. x y
##EQU00023##
[0166] The goal of the present embodiments is to minimize the
functional s[u]. Since .eta. is different from zero everywhere
except the boundaries, it is not a factor in the minimization.
Hence:
s [ u ] u = - div ( .PHI. ' ( .gradient. u ) .gradient. u
.gradient. u ) - .mu. ( f - u ) ##EQU00024##
[0167] The minimizer of s[u] satisfies the Euler-Lagrange
equation:
div ( .PHI. ' ( .gradient. u ) .gradient. u .gradient. u ) + .mu. (
f - u ) = 0 ##EQU00025##
[0168] This can be regarded as the steady-state (t.fwdarw..infin.)
of the diffusion process given by gradient descend:
u t = - s [ u ] u = div ( .PHI. ' ( .gradient. u ) .gradient. u
.gradient. u ) + .mu. ( f - u ) ##EQU00026##
where
.PHI. ' ( .gradient. u ) .gradient. u ##EQU00027##
represents a diffusion coefficient in anisotropic diffusion. the
derivative of .PHI. is:
.PHI. ' ( .gradient. u ) = .gradient. u .beta. 2 + .gradient. u 2 +
.gradient. u ##EQU00028##
[0169] Thus, the accepted total anisotropic diffusive process
is:
u t = div ( ( 1 .beta. 2 + .gradient. u 2 + ) .gradient. u ) + .mu.
( f - u ) = ( div .gradient. u .beta. 2 + .gradient. u 2 ) +
.gradient. 2 u + .mu. ( f - u ) ##EQU00029##
[0170] The chosen boundary condition for the above process is
Newman condition (mirror reflecting condition):
u.sub.n|.sub..differential..OMEGA.=0
Discretization of the Process
[0171] Following is a description for calculating the first summand
of the above equation.
Denote c ( .gradient. u ) = 1 .beta. 2 + .gradient. u 2 . div ( c (
.gradient. u ) .gradient. u ) = x ( c ( .gradient. u ) u x ) + y (
c ( .gradient. u ) u y ) = = x ( c ( .gradient. u ) ) u x + c (
.gradient. u ) u xx + y ( c ( .gradient. u ) ) u y + c ( .gradient.
u ) u yy ##EQU00030##
[0172] Calculation of the first summand of the above equation:
x ( c ( .gradient. u ) ) = c ( .gradient. u ) .gradient. u
.gradient. u x = = c ( .gradient. u ) .gradient. u u x 2 + u y 2 x
= = c ( .gradient. u ) .gradient. u ( u x 2 + u y 2 u x u xx + u x
2 + u y 2 u y u yx ) = = .gradient. u ( .beta. 2 + .gradient. u 2 )
3 / 2 ( u x u xx u x 2 + u y 2 + u y u yx u x 2 + u y 2 ) = = - u x
u xx + u y u yx ( .beta. 2 + .gradient. u 2 ) 3 / 2
##EQU00031##
[0173] Calculation of the third summand of the above equation:
y ( c ( .gradient. u ) ) = c ( .gradient. u ) .gradient. u
.gradient. u y = = c ( .gradient. u ) .gradient. u u x 2 + u y 2 y
= = c ( .gradient. u ) .gradient. u ( u x 2 + u y 2 u y u yy + u x
2 + u y 2 u x u yx ) = = - .gradient. u ( .beta. 2 + .gradient. u 2
) 3 / 2 ( u y u yy u x 2 + u y 2 + u x u yx u x 2 + u y 2 ) = = - u
y u yy + u x u yx ( .beta. 2 + .gradient. u 2 ) 3 / 2
##EQU00032##
[0174] Using the above div(c(|.gradient.u|).gradient.u)
becomes:
div ( c ( .gradient. u ) .gradient. u ) = u x u xx + u y u yx (
.beta. 2 + .gradient. u 2 ) 3 / 2 u x + c ( .gradient. u ) u xx - u
y u yy + u x u yx ( .beta. 2 + .gradient. u 2 ) 3 / 2 u y + c (
.gradient. u ) u yy = = - u x 2 u xx + u x u y u yx ( .beta. 2 +
.gradient. u 2 ) 3 / 2 + u xx .beta. 2 + .gradient. u 2 - u y 2 u
yy + u x u y u yx ( .beta. 2 + .gradient. u 2 ) 3 / 2 + u yy .beta.
2 + .gradient. u 2 = - u x 2 u xx - 2 u x u y u yx - u y 2 u yy + u
xx ( .beta. 2 + .gradient. u 2 ) + u yy ( .beta. 2 + .gradient. u 2
) ( .beta. 2 + .gradient. u 2 ) 3 / 2 = u xx ( .beta. 2 +
.gradient. u 2 - u x 2 ) - 2 u x u y u yx + u yy ( .beta. 2 +
.gradient. u 2 - u y 2 ) ( .beta. 2 + .gradient. u 2 ) 3 / 2 = = u
xx ( .beta. 2 + u y 2 ) - 2 u x u y u yx + u yy ( .beta. 2 + u x 2
) ( .beta. 2 + u x 2 + u y 2 ) 3 / 2 ##EQU00033##
[0175] The total discrete diffusive process is:
u n + 1 = u n + .DELTA. t ( u xx n ( .beta. 2 + u y n 2 ) - 2 u x n
u y n u yx n + u yy n ( .beta. 2 + u x n 2 ) ( .beta. 2 + u x n 2 +
u y n 2 ) + ( u xx n + u yy n ) + .mu. ( f - u n ) )
##EQU00034##
The first derivatives u.sub.x and u.sub.y can be calculated by
central derivative approximations:
u x , i j = u i , j + 1 - u i , j - 1 2 dx = 1 2 ( u i , j + 1 - u
i , j - 1 ) ##EQU00035## u y , i j = u i + 1 , j - u i - 1 , j 2 dy
= 1 2 ( u i + 1 , j - u i - 1 , j ) ##EQU00035.2##
[0176] The second derivatives u.sub.xx, u.sub.yy and u.sub.xy can
be calculated as follows:
u.sub.xx,ij=u.sub.i,j+1-2u.sub.i,j+u.sub.i,j-1
u.sub.yy,ij=u.sub.i+1,j-2u.sub.i,j+u.sub.i-1,j
u.sub.xx,ij=1/4(u.sub.i-1,j-1-u.sub.i-1,j+1-u.sub.i+1,j-1+u.sub.i+1,j+1)
[0177] The parameter .mu. is calculated iteratively from
Euler-Lagrange equation:
div ( .PHI. ' ( .gradient. u ) .gradient. u .gradient. u ) + .mu. (
f - u ) = 0 ##EQU00036##
[0178] The equation is multiplied by (f-u) and integrated over the
image area:
.intg. .OMEGA. ( f - u ) div ( .PHI. ' ( .gradient. u ) .gradient.
u .gradient. u ) x y = - .mu. .intg. .OMEGA. ( f - u ) 2 x y = -
.mu..sigma. 2 .intg. .OMEGA. x y ##EQU00037## .mu. = 1 .sigma. 2
.intg. .OMEGA. x y .intg. .OMEGA. ( u - f ) div ( .PHI. ' (
.gradient. u ) .gradient. u .gradient. u ) x y ##EQU00037.2##
[0179] The discrete form of .mu. is:
.mu. n = 1 .sigma. 2 M N i j ( .mu. n - f ) ( u xx n ( .beta. 2 + u
y n 2 ) - 2 u x n u y n u yx n + u yy n ( .beta. 2 + u x n 2 ) (
.beta. 2 + u x n 2 + u y n 2 ) 3 / 2 + ( u xx n + u yy n ) )
##EQU00038##
where M and N are the sizes of the image.
[0180] Although the invention has been described in conjunction
with specific embodiments thereof, it is evident that many
alternatives, modifications and variations will be apparent to
those skilled in the art. Accordingly, it is intended to embrace
all such alternatives, modifications and variations that fall
within the spirit and broad scope of the appended claims.
[0181] All publications, patents and patent applications mentioned
in this specification are herein incorporated in their entirety by
reference into the specification, to the same extent as if each
individual publication, patent or patent application was
specifically and individually indicated to be incorporated herein
by reference. In addition, citation or identification of any
reference in this application shall not be construed as an
admission that such reference is available as prior art to the
present invention. To the extent that section headings are used,
they should not be construed as necessarily limiting.
* * * * *