FIELD OF THE INVENTION

[0001] The present disclosure relates generally to image processing. More particularly, the present disclosure relates to image processing related to medical images.

BACKGROUND OF THE INVENTION

[0002] Medical imaging is used as a diagnostic tool as well as an experimental tool to study the anatomy and physiology of humans and other animals. It is also used to provide targeted treatment for some illnesses, for example, cancer. Various medical imaging techniques include X-ray; ultrasound; radiation therapy; positron emission tomography (PET); magnetic resonance imaging (MRI); and computed tomography (CT).

[0003] Digital images obtained through various medical imaging techniques are processed to obtain anatomical and physiological information. One such process is known as segmentation, in which the digital image is partitioned into multiple segments to locate various objects and boundaries of interest within the image (so that a region of interest may be extracted). Typically, segmentation involves assigning a label to every pixel in an image such that pixels with the same label share certain visual characteristics, which can then be used to identify various regions within the image having similar characteristics or computed property such as color, intensity, or texture.

[0004] Accurate segmentation of medical images is required for 3D visualization and quantification of anatomic structures. To segment an object, a method identifies all pixels (a single point in a 2D plane) or voxels (a single point in a 3D space) that bound or belong to a region of interest. The method classifies pixels or voxels based on similar characteristics such as pixel intensities, image gradients or textures.

[0005] Thresholding is a classic method for segmenting images. There are two types of thresholding methods: global and local. Global thresholding uses a constant threshold value for the entire image whereas local thresholding uses different threshold values for different regions (or objects) defined in an image. Both types of thresholding can use either single-level or multi-level thresholding depending on the pixel intensity value of the region of interest.

[0006] Single-level thresholding uses a threshold value T to turn an image l(x, y) into a binary image B(x, y) by using the equation, Γΐ, for I(x, y) > T

0, forI(x, y) >= T

[0007] Eq. 1 segments the region of interest by rejecting all pixels with intensity values lower than the threshold value T, while retaining those with intensities values higher than the threshold.

[0008] Multi-level thresholding uses two or more threshold values T1 and T2 to turn an image l(x, y) into a binary image B(x, y) by using the equation,

il, for TK I^ y) < T2

B(x,y) -

[0, otherwise

[0009] Eq. 2 segments the region of interest by rejecting all pixels with intensity values lower than the threshold value T1 and higher than the threshold value 12, while retaining those with intensities values between threshold values T1 and T2.

[00010] The application of thresholding requires knowledge about the pixel intensity values in the region of interest. The success of this method depends on choosing an appropriate threshold value. Current state-of-the-art algorithms in the field of thresholding are limited to either using fixed thresholds or analyzing the histogram at a global scale (Otsu's method), (N. Otsu, "A threshold selection method from gray-level histograms," IEEE Trans, on System, Man and Cybernetics, vol. SMC-9, pp. 62-66, 1979).

[00011] Fixed thresholds are simple to implement, but rarely work in medical images due to varying contrast across different imaging scanners and acquisition techniques used in different clinics. Complex anatomic structures and the high noise nature of medical images further compound the challenge.

[00012] Otsu's method automates the selection of thresholds by using the probability distribution function of the histogram, which shows the number of times that a pixel intensity value appears in an image. While Otsu's method works well for images with a distinct foreground object and its background, in practice, the separation of the foreground object and the background may not be clearly defined and will generally overlap each other.

[00013] Furthermore, Otsu's method has a major drawback of being able to threshold only a single foreground object; in typical medical images, many foreground objects (anatomic structures) exist. Therefore, without prior knowledge of the histogram distribution of the individual objects, appropriate threshold values are difficult to find, even by using Otsu's method. [00014] Thus, there is still a need for a method of automatically identifying complex anatomical structures in a medical image.

SUMMARY OF THE INVENTION

[00015] Herein is disclosed a histogram analysis method that combines both knowledge based intensity information and recursive subdivision (multi-scale analysis) to adaptively find thresholds in medical images. The medical images may be from different scanners and/or clinics.

[00016] In a first aspect, the present disclosure provides a histogram analysis method for automatically identifying a desired anatomic structure in a medical image. The method comprises dividing a histogram associated with the medical image to form at least one histogram segment having the desired anatomical structure based on expected intensity values of the desired anatomic structure; detecting one or more peak intensity values within the histogram segment; recursively sub-dividing the histogram segment upon detecting more than one peak intensity value within the histogram segment; comparing each histogram segment with a peak intensity value to the expected intensity values of the desired anatomical structure; and automatically identifying the histogram segment with the intensity value associated with the desired anatomic structure in the histogram segment containing only one peak intensity value.

[00017] In a further aspect the method comprises segmenting the anatomic structure by thresholding the identified intensity values of the histogram segment associated with the anatomic structure to form a segmented image of the desired anatomic structure.

[00018] In a further aspect, the method comprises detecting one or more expected intensity values in the array of integers corresponding to the histogram segment for the desired anatomic structure; recursively sub-dividing the array of integers and recursively sub-dividing the histogram upon detecting more than one expected intensity value in the array of integers.

[00019] In a further aspect, the step of sub-dividing the array of integers and the histogram comprises: computing a maximum separation range between the more than one expected intensity values; cropping the histogram in the determined maximum separation range; determining a deepest trough in the cropped histogram; and subdividing the array of integers and the histogram based on the deepest trough.

[00020] In a further aspect, the present disclosure provides a system for automatically identifying a desired anatomic structure in a medical image, the system comprising a processor adapted to perform a method for automatically identifying an anatomic structure as disclosed herein.

[00021] In a further aspect, the present disclosure provides a computer program product, readable by a computer and containing instructions operable by a processor of a computer system to cause the processor to perform a method as disclosed herein.

[00022] Other aspects and features of the present disclosure will become apparent to those ordinarily skilled in the art upon review of the following description of example embodiments in conjunction with the accompanying figures. BRIEF DESCRIPTION OF THE DRAWINGS

[00023] Example embodiments will now be described, by way of example only, with reference to the attached Figures, wherein:

[00024] Fig. 1 is a flowchart of a histogram analysis method in accordance with an aspect of the present disclosure;

[00025] Fig. 2 is a flowchart of a histogram analysis method in accordance with another aspect of the present disclosure;

[00026] Fig. 3 is a flowchart of a peak detection technique.

[00027] Fig. 4 is an illustration of histograms of images obtained from three different scanners analyzed by a histogram analysis method of the present disclosure; and

[00028] Fig. 5 is a figure showing heart segmentation of computed tomography angiography (CTA) images obtained from three different scanners using a histogram analysis method of the present disclosure. DETAILED DESCRIPTION

[00029] In accordance with the present disclosure, a knowledge-based multi-scale histogram analysis method is provided that combines: (1 ) multi-scale histogram analysis and (2) a priori knowledge to adaptively calculate the appropriate thresholds for anatomic structures in a broad range of medical images.

[00030] Generally, the present invention provides a method and system for automatically identifying a desired anatomic structure in a medical image, for example, within a computed tomography (CT) image. For ease of understanding, the method and system for automatically identifying a desired anatomic structure in a medical image will be described with the example of cardiac CT images. However, the methods described herein are not limited to CT images and may be applicable to medical images in general. For example, the image may be from an X-ray; ultrasound; radiation therapy; positron emission tomography (PET); or magnetic resonance imaging (MRI).

[00031] As used herein, the term "recursive" refers to a sub-dividing step that can repeat itself indefinitely. The recursive method calls itself to decide if further detection and division is needed. "Recursively" is used interchangeably with "iterative". The term "iterative" is used herein to mean repeating any number of times as required.

[00032] As used herein, the term "expected intensity values" refers to the intensity value or values of pixels in a knowledge base associated with a particular anatomic structure.

[00033] As used herein, the term "array of integers" refers to a priori information of the expected pixel intensity values for different anatomic structures in an image or volume. The array of integers is stored in a knowledge base. The term knowledge base is used interchangeably with the term database.

[00034] A knowledge base may be provided in the form of a user input or computed from a statistical background. It contains existing information on the expected peak intensity values for anatomic structures. As can be seen in Fig. 4, different types of scanners show slightly different peaks for the same anatomical structure. This a priori knowledge is stored in the knowledge base and is used to drive the subdivision of the histogram and appropriately label anatomical structures from different image sources. The knowledge base may contain any intensity values or values suitable for the analysis of digital images. For example, when analyzing CT images a knowledge base containing intensity values expressed in Hounsfield Units may be used.

[00035] Fig. 1 shows a method for automatically identifying a desired anatomic structure in a medical image (100) according to an aspect of the present disclosure. At 102, a histogram associated with the medical image is divided to remove unwanted anatomic structures and to form at least one segment having the desired anatomical structure based on an expected intensity value of the desired anatomic structure. At 104, for each histogram segment, one or more peak intensity value is detected. If more than one peak intensity value is detected, the histogram segment is recursively sub-divided into further histogram segments. If only one peak intensity value is detected, at 106, the peak intensity value of the histogram segment is compared to the expected intensity value of the desired anatomic structure. At 108, the anatomic structure associated with the peak intensity value in the histogram segment containing only one peak intensity value is automatically identified. [00036] Example 1

[00037] Given an image or a volume of images, the knowledge-based multi-scale histogram analysis method iteratively divides the histogram into two histograms using a priori information provided by the knowledge base. The histogram may be split into more than two histograms. The value of two is used solely as an example. A flowchart of a method of histogram analysis for automatically identifying a desired anatomic structure is shown at Figure 2.

[00038] In the knowledge base used in this example, lungs have a Hounsfield Unit (HU) of approximately 800 while muscles have a HU of approximately 100. In the pseudocode shown in Table 1 , the knowledge base is provided as a Function 1 (MultiscaleHistogramAnalysis) via the input variable knowledge, which is an array of integers containing the expected pixel intensities for different anatomic structures. Table 1 outlines Function 1. Table 1

Function 1 : Performs knowledge-based multi-scale histogram analysis

function MultiscaleHistogramAnalysis(histogram, knowledge, output)

{

// Remove histogram values from unwanted anatomic structures

sort knowledge in ascending ordered

set minValue to percentA of the minimum value in knowledge

set maxValue to percentB of the maximum value in knowledge

for each item_i in histogram

if i < minValue or i > maxValue

set itemj to 0

// Perform recursive multi-scale histogram analysis

AnalyzeHistogramSegment (histogram, knowledge, output)

}

[00039] At 202, medical images are obtained and a histogram is computed at 204. Optionally, at action 206, the method performs preprocessing of the histogram by removing pixel values from unwanted anatomic structures, for example, as specified in Function 1 of the pseudocode shown in Table 1. Unwanted anatomic structures are determined as any pixel intensities falling outside a specified range of pixel intensities provided by the knowledge base. Any pixel intensities that are not within the range are removed leaving only the desired pixel intensities. The cutoff variables used in Table 1 were, for example, 90% for percentA and 110% for percentB. However any suitable cutoff range may be used, depending on the type of images being analyzed. The step of preprocessing is optional and may occur throughout the method.

[00040] Segmentation

[00041] At 208 the method divides the histogram into histogram segments based on the pixel intensity of anatomic structures provided in the knowledge base. Given a set of pixel intensity values in the knowledge base, the method first separates the values into two groups by finding the maximum spread between two adjacent values. The mean pixel intensity is the calculated for each group. Using the mean two pixel intensity values, la and lb (corresponding to each group), the method finds the deepest trough in the histogram between the range [la, lb] as the dividing point between the two histograms.

[00042] At 210, the method analyzes the knowledge base and determines the number of values in the knowledge base at 212. If only one value is found, the method performs a peak detection technique to find the peak that is closest to the one specified in the knowledge base at 214. This peak is labeled as the anatomic structure specified in the knowledge base at 216, by assigning the output variable with the intensity value of the peak, and the segmentation analysis is complete.

[00043] If more than one value is found in knowledge, the method splits the histogram and knowledge into two segments, histograml , histogram2, knowledge^ knowledge2 at 224, by first computing the maximum separation between the values in knowledge at 218, followed by detecting the deepest trough (the lowest point or valley) in the maximum separation range at 220. The trough is used to split both histogram and knowledge at 222. This process continues until only one value is found in the sub-divided knowledge segment. For each of these segments, the method then performs the same peak detection and peak labeling technique as described in the previous paragraph, which produces an output variable containing the intensity value of the identified anatomic structure. The resulting output variables from each histogram segment (due to subdivision) are then concatenated into a single output array. This produces a final output array containing the identified pixel intensities for each desired anatomic structure in the knowledge base. The method continues iteratively until all the peaks are labeled.

[00044] Table 2 outlines the recursive histogram analysis performed in Function 2 (AnalyzeHistogramSegment). Table 2

Function 2: Performs recursive histogram segment analysis

function AnalyzeHistogramSegment(histogram, knowledge, output)

{

if the number of items in knowledge is equal to one

DetectPeaks(histogram, peakValues, peakLocations, troughValues, troughLocations)

// Get the peak that is closest to the item in knowledge

for each itemj from 1 to N in peakValues

compute distances,

where distancej = item in knowledge - itemj in peakValues set output to the minimum of all distances

return

else

// Compute the maximum separation between values in knowledge set index to 0

set maxLocation to 0

set maxSeparation to -infinity

for each item j from 1 to N in knowledge

compute separationj = item_i - the next item in knowledge, item_(i+1 )

if separationj is greater than maxSeparation

set maxLocation to index

set maxSeparation to separationj

set index to index + 1

// Crop the histogram in the area of maximum separation

set startSe pa ration to mean of item l to item maxLocation in knowledge set endSeparation to mean of item_(maxLocation+1 )to itemjast in knowledge

set newHistogram to the segment between

startSeparation and endSeparation in histogram

// Detect the deepest trough in the separation range

DetectPeaks(newHistogram, peakValues, peakLocations, troughValues, troughLocations)

set index to 0

set minPosition to 0

set minValue to -infinity

for each item j from 1 to N in troughValues

if itemj is less than minPosition

set minValue to itemj

set minPosition to index

set index to index + 1

compute splitlndex = minPosition + startSeparation // Split histogram and knowledge and perform recursive analysis by a sub- scale

Split histogram into histograml and histogram2,

where histograml is the segment from 0 to splitlndex and histogram2 is the segment from splitlndex + 1 to the end of histogram

Split knowlege into knowledgel and knowledge2,

where knowledgel is the segment from 0 to maxLocation and knowledge2 is the segment from maxLocation + 1 to the end of knowledge

AnalyzeHistogramSegment (histograml , knowledgel , outputl )

AnalyzeHistogramSegment (histogram2, knowledge2, output2) concatenante outputl and output2 into output

return

}

[00045] Fig. 3 outlines the steps that are involved in the peak detection technique.

Function 3 detects peaks and troughs in a histogram by tracking the changes in the values of the histogram. Given a histogram, histogram = {hi , h2... hn}, in the form of an array of values, a peak in the histogram is defined as a point, hp, whenever there exists two other points, ha, hb, with the conditions of a < p < b and ha, hb < hp. Similarly, a trough in the histogram is defined as a point, ht, whenever there exists two other points, he, hd, with the conditions of c < t < d and he, hd > ht.

[00046] To detect all the peaks and troughs as defined by the previous two definitions, the peak detection technique, for example, steps through all the values in the histogram in the order of 1 to N at 302. At each value of the histogram, hi, the technique keeps track of the current maximum value, currentMax and its corresponding location, maxPosition, as well as the current minimum value, currentMin and the its corresponding location, minPosition, by comparing the current value, hi, to the currentMax and currentMin. If the current value is greater than the current maximum, then the current maximum becomes the current value at 304. Similarly, if the current value is smaller than the current minimum, then the current minimum becomes the current value at 306.

[00047] If the current value, hi, is smaller than the current maximum, currentMax, by a preset amount of peakDeltal , then a peak is detected at 308. The location and value of the peak is added to the outputs of peakValues and peakLocations. If the current value, hi, is greater than the current minimum, currentMin, by a preset amount of peakDelta2, then a trough is detected at 310. The location and value of the trough is added to the outputs of troughValues and troughLocations. The selection of peakDeltal and peakDelta2 can be varied based on the type of histogram being analyzed. In the example illustrated in Fig. 5, the peakDeltal and peakDelta2 is set to 10% of the maximum value of the histogram.

[00048] The steps involved in the detection of the location and values of all peaks and troughs in the histogram using Function 3 are outlined in Table 3.

Table 3

Function 3; Detect the locations and values of ail peaks and troughs in the hislogram function DetectPeaks(histogram, peakValues, peakLocations, troughValues, troughLocations)

{

set index to 0

set minPosition to 0

set maxPosition to 0

set currentMax to -infinity

set currentMin to infinity

set findMax to true for each item from 1 to N in histogram

if item_i is larger than currentMax

set currentMax to item J

set maxPosition to index

if item i is smaller than currentMin

set currentMin to item J

set minPosition to index

// If in the mode of looking for peaks

if findMax is true

if item_i is less than currentMax - peakDelta

add item_i to peakValues

add index to peakLocations

set currentMin to itemj

set minPosition to index

set findMax to false

// Else in the mode of looking for troughs

else if findMax is false

if itemj is greater than currentMin + peakDelta

add item_i to troughValues

add index to troughLocations

set currentMax to itemj

set maxPosition to index

set findMax to true

set index to index + 1

} [00049] Fig. 4 shows the histogram analyses of three Computed Tomography Angiography (CTA) images generated by three different scanners (GE, Phillips and Siemens). The x axes of the graphs show pixel intensities, while the y axes show pixel counts. As can be seen in Fig. 4 these histograms are very dissimilar. Using a histogram analysis method described herein the histogram peaks are appropriately labeled as lungs, muscles and vessels by computing intensity thresholds.

[00050] Once an intensity value is automatically identified, a particular anatomic structure may be selected from the histogram using automatic thresholding. The automatic thresholding using the intensity value will extract the anatomic structure from the image/volume. The method of thresholding removes all intensity values not associated with the anatomic structure in the image/volume, leaving only the desired intensity values that are associated with the anatomic structure.

[00051] Fig. 5 shows thresholding using a histogram analysis method as described above. In this example the knowledge base used is a knowledge base of Hounsfield units (HU) as follows:

Lungs: ~ -800 HU

Fat: ~ -200 HU

Muscles: ~ +100 HU

Vessels: - + 400 HU

[00052] The method is able to automatically compute the appropriate thresholds for various anatomic structures. In the top row, original CT images from three different imaging systems are shown. Segmentation of the original images using the method is shown in the middle row. The bottom row shows a three dimensional reconstruction of the segmented CT volumes. Although the original images (top row) have vastly different contrast levels, the method was able to correctly find the appropriate threshold for the heart and segment the heart correctly.

[00053] An embodiment is a system for automatically identifying an anatomic structure in a medical image. The system comprises a processor adapted to perform a histogram analysis method as described herein.

[00054] An embodiment is a computer program product, readable by a computer and containing instructions operable by a processor of a computer system to cause the processor to perform a histogram analysis method for automatically identifying a desired anatomic structure in a medical image.

[00055] The computer program product can be represented as a software product stored in a machine-readable medium (also referred to as a computer-readable medium, a processor-readable medium, or a computer usable medium having a computer- readable program code embodied therein). The machine-readable medium can be any suitable tangible medium, including magnetic, optical, or electrical storage medium including a diskette, compact disk read only memory (CD-ROM), memory device (volatile or non-volatile), or similar storage mechanism. The machine-readable medium can contain various sets of instructions, code sequences, configuration information, or other data, which, when executed, cause a processor to perform steps in a method according to an embodiment of the invention. Those of ordinary skill in the art will appreciate that other instructions and operations necessary to implement the described invention can also be stored on the machine-readable medium. Software running from the machine- readable medium can interface with circuitry to perform the described tasks.

[00056] In the preceding description, for purposes of explanation, numerous details are set forth in order to provide a thorough understanding of the embodiments of the invention. However, it will be apparent to one skilled in the art that these specific details are not required in order to practice the invention. In other instances, well-known electrical structures and circuits are shown in block diagram form in order not to obscure the invention. For example, specific details are not provided as to whether the embodiments of the invention described herein are implemented as a software routine, hardware circuit, firmware, or a combination thereof.

[00057] The above-described embodiments of the invention are intended to be examples only. Alterations, modifications and variations can be effected to the particular embodiments by those of skill in the art without departing from the scope of the invention, which is defined solely by the claims appended hereto.