Automated macromolecular crystal detection system and method

[0001] This application claims priority in provisional application filed on May 30, 2002, entitled “Augmented Automated Macromolecular Crystal Detection from Light Microscopy Images” Ser. No. 60/385,210, by inventors Christian et al.

[0002] The United States Government has rights in this invention pursuant to Contract No. W-7405-ENG-48 between the United States Department of Energy and the University of California for the operation of Lawrence Livermore National Laboratory.

[0003] The present invention relates to automated image recognition and detection systems and methods, and more particularly to an automated system and method using a phase-based edge detection process to detect macromolecular crystals from two-dimensional images obtained from light microscopy of crystallization experiments.

[0004] Proteomics is the field of bioscience involving the characterization of the proteins encoded by the human genome, and enabled by the gene sequence data produced by the human genome project. Since the structure of a protein is key to understanding its function, one field of proteomics in particular has rapidly emerged concerning high throughput structure determination or structural genomics. In determining protein structure, the proteins are first crystallized, and then an X-ray generator produces diffraction patterns from which a three-dimensional picture of the atomic arrangement in the crystal can be obtained. Advances in macromolecular crystallography techniques, computer speed, and the availability of high-energy synchrotron x-ray sources, make rapid structure determination possible given adequate quality protein crystals.

[0005] Crystal growth, however, is difficult because proteins are large, irregularly shaped molecules that do not readily come together in a repeating pattern, and the complete set of crystallization conditions is too large and impractical to screen comprehensively. Thus, previously uncrystallized proteins must be screened on a trial and error basis against a large array of conditions that have the potential to induce crystal formation. Automated methods using, for example, robotic liquid handling devices, robotic CCD-based microscope cameras, or light microscopes equipped with robotic stages and CCD cameras, have been developed and are commercially available to speed up the process of setting up and recording the results (automated image capture) of a large number of crystallization trials. However, a practical problem remains in that each experiment must still be visually inspected to determine successful crystal formation. In fact, the high throughput enabled by the automation in setup and image-capture has increased the visual inspection bottleneck, which is typically performed manually by human intervention.

[0006] One example of an automated crystal detection method developed to address the visual inspection bottleneck is disclosed in the article “Intelligent Decision Support for Protein Crystal Growth” (by Jurisica et al, IBM Systems Journal, Vol. 40, No. 2, 2001). In that article, and as shown in FIGS. 3–7 thereof, images of screening results are analyzed using a two-dimensional Fourier transform. In particular, FIG. 5C illustrates the Fourier frequency spectrum used in the analysis, and FIG. 5D illustrates an analysis of the spectrum derivatives and circular averages to provide features information of the image. From this feature extraction and analysis, the outcome of the experiment is classified as, clear drop, amorphous precipitate, phase separation, microcrystals, crystals, or unknown.

[0007] Despite such efforts, difficulties in automating (i.e. without human intervention) crystal detection remain due to such factors as poor image quality due to noise and low contrast, differences in crystal shapes, poorly formed crystals, etc. With respect to poor image quality, crystals may have less contrast relative to the background than other objects or particles. For example, the difference between the crystal and the background based on 256 gray levels is often 15 levels, whereas the difference for dirt is usually above 40 levels. Additionally, many different crystal shapes exist due to, for example, the existence of several large classes of crystal shape, the picture is a 2-D projection of a 3-D object, crystal imperfections with faulty edges, and large variations in crystal size, e.g. ranging from about 10 μm to nearly 300 μm. There are also many things on the picture that are not real crystals, such as dirt, precipitation, quasi-crystals, small drop due to condensation, and unidentified effects. Additionally, an automated crystal detection process must also achieve a high threshold of accuracy by being able to identify virtually all crystals with a low false-positive rate.

[0008] Thus in summary, there is a need for an automated crystal detection method and system for inspecting two-dimensional images and successfully detecting crystals therefrom. An automated solution for crystal detection, such as implemented by a software program, would be a great labor savor by possessing the capability of processing thousands of images a day and provide analysis substantially free from false positives.

[0009] One aspect of the present invention includes a method of detecting macromolecular crystals in light microscopy images comprising: detecting edges in said images by identifying local maxima of a phase congruency-related function associated with each image; segmenting the detected edges into discrete line segments; evaluating the geometric relationships that the line segments have with each other to identify any crystal-like qualities; and determining the presence of crystals in each image based on said evaluation.

[0010] Another aspect of the present invention includes a computerized system for detecting macromolecular crystals from light microscopy images comprising: a digital conversion component that converts said light microscopy images into corresponding phase-based digital image data using the Fourier transform; an edge detection component that detects edges from the image data by computing local maxima of a phase congruency-related function associated with each image; a segmentation component that divides the detected edges into discrete line segments; and a geometric analyzer component that evaluates the geometric relationships that the line segments have with each other to identify any crystal-like qualities, and determines whether crystals are present in each image based on said evaluation.

[0011] Another aspect of the present invention includes a computerized system for detecting macromolecular crystals from light microscopy images comprising: means for digitally converting said light microscopy images into corresponding phase-based digital image data using the Fourier transform; means for detecting edges from the image data by computing local maxima of a phase congruency-related function associated with each image; means for dividing the detected edges into discrete line segments; means for evaluating the geometric relationships that the line segments have with each other to identify any crystal-like qualities; and means for determining the presence of crystals in an image from said evaluation.

[0012] Another aspect of the present invention includes a computer program product comprising: a computer useable medium having a computer readable code embodied therein for causing the detection of macromolecular crystals in light microscopy images, said computer program product having: computer readable program code means for causing a computer to detect edges in said images by identifying local maxima of a phase congruency-related function associated with each image; computer readable program code means for causing said computer to segment the detected edges into discrete line segments; computer readable program code means for causing said computer to evaluate the geometric relationships that the line segments have with each other to identify any crystal-like qualities; and computer readable program code means for causing said computer to determine the presence of crystals in each image based on said evaluation.

[0013] Another aspect of the present invention includes an article of manufacture comprising: a computer useable medium having a computer readable code means embodied therein for causing the detection of macromolecular crystals in light microscopy images, said computer readable code means in said article of manufacture comprising: computer readable program code means for causing a computer to detect edges in said images by identifying local maxima of a phase congruency-related function associated with each image; computer readable program code means for causing said computer to segment the detected edges into discrete line segments; computer readable program code means for causing said computer to evaluate the geometric relationships that the line segments have with each other to identify any crystal-like qualities; and computer readable program code means for causing said computer to determine the presence of crystals in each image based on said evaluation.

[0014] And another aspect of the present invention includes a program storage device readable by a machine, tangibly embodying a program of instructions executable by the machine to perform method steps for detecting macromolecular crystals from light microscopy images, said method steps comprising: detecting edges in said images by identifying local maxima of a phase congruency-related function associated with each image; segmenting the detected edges into discrete line segments; evaluating the geometric relationships that the line segments have with each other to identify any crystal-like qualities; and determining the presence of crystals in each image based on said evaluation.

[0024] The present invention is an automated method and system (hereinafter “method”) for detecting the presence of macromolecular crystals from light microscopy images, such as those obtained from crystallography experiments. Generally, the automated method utilizes the phase information of the pixels in each image in performing crystal detection. This provides greater resistance to noise and other artifacts than amplitude based properties, e.g. pixel intensity, which can lead to poor image quality. And in particular, the phase information of the pixels is used to identify specific geometric features in an image attributable to and most likely indicative of a crystalline structure, such as, for example, parallel edges, of similar length, facing each other, in relatively close proximity. Evaluation of detected geometries in this manner makes possible a low rate of false-positive detections and an effective resolution to the visual inspection bottleneck discussed in the Background caused by manual inspections. It is notable that the term “automated” suggests the absence of human intervention, although oversight of detection performance may still be present.

[0025] Implementation of the present invention is suitably achieved using, for example, software, computer code, ROM, integrated circuit, etc. (hereinafter “software”) to execute and control the method steps and system functions. It is appreciated that the software may be written in any suitable programming language for operation on a suitable operating system or platform, not limited to any particular language or operating system or platform. For example, the software may be written as object oriented code having implementable subroutines that enable the user to call them in sequence to solve problems.

[0026] Turning now to the figures, FIG. 1 shows an overview flow diagram of the steps involved in a first general embodiment of the present invention, indicated at reference character 10. Two-dimensional images captured from crystallography experiments are first provided as input 11 to begin the process. At reference character 12, edge features, if any, in each image are detected using phase congruency, and in particular, a phase congruency-based function. The detected edges are then segmented into discrete, straight line segments in the line segmentation step of 13, followed by a geometric analysis and evaluation (hereinafter “evaluation”) of the geometric relationships that the line segments have with each other to identify any crystal-like qualities, at 14. In particular, geometric relationships, such as parallel lines, facing each other, and similarity in length, as well as relative proximity to each other, are all crystal-like qualities in that they are characteristic of crystalline structure. And finally at 15, a determination and decision is made based on the preceding geometric evaluation, as to whether crystals are present in a particular image. A decision with respect to any one image may be reported separately, or together in a list with other images determined to have detected crystals. The corresponding ones of the experiments may subsequently be collected for X-ray determination of atomic structure and protein identification.

[0027] Generally, the edge detection stage 12 uses the phase of the Fourier transform of the image to find the images. Taking, for example, a 1-D step signal (representing the equivalent of an edge) and its cosine decomposition shown in FIG. 5, the phases of all cosines are equal to zero only at the edge location, while everywhere else, the phases have different values. Since the phases of real signals are never actually equal, phase congruency or repartition is a method for measuring how much the phases of all the cosines are equal, with an edge corresponding to a local maximum of the phase congruency. It is preferred to compute the standard deviation of the cosine of the phases since computing the standard deviation of the phases at each point can give wrong results for the simple reason that 0=2π. The phase congruency function for a 1-D continuous signal is therefore defined as:

[0028] $\begin{array}{cc}\mathrm{PC}\left(x\right)=\underset{\theta \in \left[0,2\pi \right]}{\max }\frac{\int a_{\omega }\cos \left(\omega x+{\Phi }_{\omega }-\theta \right)d\omega }{\int a_{\omega }d\omega }& \left(1\right)\end{array}$
where α_ω and Φ_ω are respectively the amplitude and phase of the signal. Thus identifying the local maxima of the phase congruency function will also identify and thereby detect the edges. It is notable, however, that while this function provides the correct computational results, it is slow to compute directly.

[0029] In the alternative, the same result can be obtained using local energy, and Gabor filters which are similar in function to those present in the human brain for image perception and vision. The local energy of a 1-D signal is defined as:
LE(x)=√{square root over (I²(x)+H²(x))}{square root over (I²(x)+H²(x))}  (2)
where I(x) is the 1-D signal, and H(x) is the Hilbert transform of I (i.e. ninety degree phase shift of I(x) in the frequency domain). Given the cosine decomposition of I(x) is:
∫α_ω cos(ωx+Φ_ω)dω  (3)
then H(x) has the decomposition:
−∫α_ω sin(ωx+Φ_ω)dω  (4)
I(x) is obtained by convolving the original signal by a filter to remove the DC component. And H(x) is obtained by filtering the previous result by the Hilbert transform of the first filter. The human visual system has neuronal structures similar to a pair of odd and even symmetric filters in quadrature. If M_e is the even filter and M_o is the odd filter, one has:

[0030] $\begin{array}{cc}{M}_e*f\left(x\right)\approx {\int }_{-\infty }^{+\infty }{a}_{\omega }\cos \left(\omega x+{\Phi }_{\omega }\right)d\omega & \left(5\right)\end{array}$
and

[0031] $\begin{array}{cc}{M}_o*f\left(x\right)\approx -{\int }_{-\infty }^{+\infty }{a}_{\omega }\sin \left(\omega x+{\Phi }_{\omega }\right)d\omega & \left(6\right)\end{array}$
The even symmetric filter is chosen so that it covers as much of the frequency spectrum as possible, and at the same time removing the D.C. term. Since it is nearly impossible to have such a perfect filter, the approximate equal sign is used in Equations (5) and (6) above. The odd filter is the π/2 phase shift even one. In the real implementation, the filters are band-pass filters and different local energies are computed for different scale. This allows for a multi-scale analysis of the image to be performed, with the possibility to look for features of different size, and for example getting rid of noise or features very small that cannot be crystal.

[0032] The phase congruency function and the local energy function are related by:

[0033] $\begin{array}{cc}\begin{array}{c}\mathrm{LE}\left(x\right)=({\int }_{-\infty }^{+\infty }{a}_{\omega }\cos \left(\omega x+{\Phi }_{\omega }\right)d\omega ,\\ -{\int }_{-\infty }^{+\infty }{a}_{\omega }\sin \left(\omega x+{\Phi }_{\omega }\right)d\omega \\ =\underset{0\in \left[0,2\pi \right]}{\max }({\int }_{-\infty }^{+\infty }{a}_{\omega }\cos \left(\omega x+{\Phi }_{\omega }\right)d\omega ,\\ -{\int }_{-\infty }^{+\infty }{a}_{\omega }\sin \left(\omega x+{\Phi }_{\omega }\right)d\omega )·\left(\cos \theta ,-\sin \theta \right)\\ =\underset{0\in \left[0,2\pi \right]}{\max }{\int }_{-\infty }^{+\infty }{a}_{\omega }(\cos \left(\omega x+{\Phi }_{\omega }\right)\cos \theta +\\ \sin \left(\omega x+{\Phi }_{\omega }\right)\sin \theta )d\omega \\ =\underset{0\in \left[0,2\pi \right]}{\max }{\int }_{-\infty }^{+\infty }{a}_{\omega }(\cos \left(\omega x+{\Phi }_{\omega }-\theta \right)d\omega \end{array}& \left(7\right)\\ \mathrm{LE}\left(x\right)=\int a_{\omega }d\omega ·\mathrm{PC}\left(x\right)& \left(8\right)\end{array}$
As used herein and in the claims, the term “phase congruency-based function” is used to describe both the phase congruency function PC(x) as well as the local energy function LE(x). A local maximum in the local energy corresponds to a local maximum in the phase congruency, and to an edge. Therefore, in order to search for local maxima in the phase congruency function, one equivalently searches for local maxima in the local energy function. These local maxima will occur at step edges of both parity (up or down), lines and bar edges, and other types of features such as the illusion patterns mentioned before. While the previous calculations for both PC(x) and I(x) were for 1-D signals, it is appreciated that a 2-D signal can also be decomposed in a series of 1-D signals by traversing and accounting for many orientations. The local energy of a point will be the maximum local energy among all the orientations.

[0034] FIG. 2 shows a detailed flow diagram, indicated at reference character 12′, of an exemplary edge detection step using a phase congruency-based function (e.g. local energy function) associated with an image. As indicated at reference character 20, and as discussed above, multiple orientations are considered. For each orientation and each scale at 21, a Gabor filter 22 is used to filter the image. The objective of this filtering step is to convert the spatial information in another space representation where the noise will be present everywhere, and the information will be more localized, so that the noise and information can be separated, and converted back to the normal space. Gabor filters are used since they closely approximate the human visual system, as previously discussed. They are the combination of the two real filters in quadrature, and are the result of a sine wave modulated by a Gaussian function as in:

[0035] $\begin{array}{cc}g\left(x,y\right)=\frac{1}{2{\mathrm{\pi \sigma }}^2}\exp \left\{-\frac{x^{\mathrm{\prime 2}}+y^{\prime 2}}{2{\sigma }^2}\right\}\cos \left(2\mathrm{\pi \omega }{x}^{\prime }\right)& \left(9\right)\end{array}$
where (x′, y′)=(x cos θ+y sin θ, −x sin θ+y cos θ). In other words, (x′, y′) is the θ rotation of (x, y). And at reference character 23, local energy is calculated from the information at the different scales.

[0036] At reference characters 24 and 25, noise variance and mean of the local energy is estimated, and a noise threshold is applied on the local energy. Considering first the noise threshold, the filter response will be decreased by the noise threshold to remove the components considered as noise. However, finding the right threshold can be difficult. The expected response of the Gabor filters to a pure noise signal must first be examined. If the noise is supposed to be Gaussian, the complex response will have a 2-D Gaussian distribution. And the magnitude will be a Rayleigh distribution as:

[0037] $\begin{array}{cc}R\left(x\right)=\frac{x}{{\sigma }_g^2}{e}^{\frac{-x^2}{2{\sigma }_g^2}}& \left(10\right)\end{array}$
where σ_g² is the variance of the 2-D Gaussian distribution. The mean and variance of the Rayleigh distribution are given by:

[0038] $\begin{array}{cc}{\mu }_r={\sigma }_g\sqrt{\frac{\pi }{2}}\mathrm{and}{\sigma }_r^2=\frac{4-\pi }{2}{\sigma }_g^2& \left(11\right)\end{array}$
The noise threshold can be set to some number of standard deviation beyond the mean of the distribution as in T=μ+kσ_r.

[0039] Additionally, the noise amplitude distribution must be determined. The median value is considered as a statistically robust estimation of the mean. The smallest scale of the Gabor filter is used because the noise should be the most present at this scale. The median of a Rayleigh distribution is the value x such that

[0040] $\begin{array}{cc}{\int }_0^x\frac{x}{{\sigma }_g^2}{e}^{\frac{-x^2}{2{\sigma }_r^2}}=\frac{1}{2}& \left(12\right)\end{array}$
that leads to
median=σ_g√{square root over (−2ln(1/2))}  (13)
The mean of the Rayleigh distribution being σ_g√{square root over (π/2))}, one gets the estimated mean of the amplitude response at smallest scale by a simple proportionality

[0041] $\begin{array}{cc}\mathrm{mean}=\frac{{\sigma }_g\sqrt{\pi /2}}{{\sigma }_g\sqrt{-2\ln \left(1/2\right)}}·\mathrm{median}& \left(14\right)\end{array}$
Assuming that noise has a uniform repartition among frequency, the noise threshold at larger scales can be deduced by the threshold at the smallest scale since the noise amplitude response will be proportional to the bandwidth, itself related to the scale. It is notable that while the noise reduction process is described as executed in the edge detection stage 12′, it is appreciated that noise reduction may be performed as a pre-processing step prior to edge detection.

[0042] At reference character 26, local energy from each orientation is summed together, and the magnitude from each orientation is summed together at 27. The summations provide the values necessary to solve Equation (8) above. In particular, by dividing local energy by the magnitude, as indicated at reference character 28, the phase congruency function is obtained. Furthermore, at reference character 29; local maxima of the phase-congruency-based function (either PC(x) or LE(x)) is identified to detect corresponding edges. In this manner an edge map may be constructed as shown in FIG. 6, for example, with local maxima occurring at step edges of both parity (up or down), lines and bar edges, and other types of features such as the illusion patterns.

[0043] FIG. 3 shows a flow diagram of a preferred geometry analysis and evaluation process, generally indicated at reference character 14′. As discussed with respect to FIG. 1, the geometric evaluation stage at 14 occurs subsequent to the segmentation 13 of the detected edges into discrete line segments, and evaluates the geometric relationships that the line segments have with each other to identify crystal-like qualities. As indicated at 30, each combination of pairs of segments are first rated based on characteristically crystal-like qualities, including, but not limited to: parallel lines, facing each other, similarity in length, and relative proximity to each other. Next at 31, the line segments are grouped based on the preceding crystal-like quality ratings. For each group at 32, each edge pixel is classified according to the local edge orientation at 33. For example, FIG. 7 shows a multiple screen shot view of the edge map of FIG. 5 separated into eight images, with each image containing features having a given orientation. For each orientation at 34, auto-correlation is calculated to enhance shapes with parallel lines, and subsequently analyzed at 35 to decide if the shape is a crystal or not. FIG. 8 illustrates two orientations before (left column) and after (right column) auto-correlation is performed, showing the enhancements. Additionally at 35 the quality of the auto-correlation is rated using the following criteria: numbers of features in the auto-correlation above specific threshold, and length and location of these features. A global “quality score” is then calculated from the auto-correlation quality ratings. At 36 a final decision is made from the global quality score for each orientation as to whether a crystal has or has not been detected. Both statistical analysis and neural networks are used for this purpose. Depending on such determination that a crystal has in fact been detected, the results are written to a crystal list at 37, which may subsequently be viewed or otherwise reported to a researcher, scientist, or operator, or stored for later retrieval. The process may also be augmented by ranking the images according to crystal quality as determined from the aforementioned steps.

[0044] Additionally, the geometry analysis and evaluation stage may also be performed using additional criteria, including: (1) minimum gradient along edges, and (2) central symmetry. With respect to the first additional criteria, and for a given edge, the gradient (i.e. first derivative of the original picture) for each pixel overlapping this edge is summed together. If this sum is above a specific gradient, the edge is then considered relevant. With respect to the second criteria, symmetry centers are found using the phase information. An important aspect of symmetry is the periodicity that it implies in the structure of the object. To determine the centers of symmetry, the picture is filtered with the same pair of even and odd symmetric Gabor filters as for the edge detection technique. At a point of symmetry, the absolute value of the even filter will be large, and the absolute value of the odd one will be small. Taking the difference of the absolute value of the even-symmetric filter and the odd one give a quantification of the degree of symmetry in an object.

[0045] A flow diagram of an exemplary second embodiment is also shown in FIG. 4. Similar to the general algorithm of FIG. 1, the process here also begins with the input of images at 41, captured from, for example, crystallography experiments. Steps 42 and 43 illustrate optional pre-processing functions available for improving the crystal detection process. At 42, noise reduction is performed as a pre-processing step, similar to that previously discussed. Additionally, at step 43, a smoothing filter is employed, such as that performing nonlinear anisotropic diffusion, to clean the image prior to edge detection. Nonlinear anisotropic diffusion filters are commonly known and used for image enhancement while preserving edges. In particular, the filter performs intraregion smoothing while inhibiting interregion smoothing. It is notable also that image cleaning by smoothing step may be alternatively employed following the edge detection process, to clean and smooth the edge map of an image (not shown). At 44 edge detection is performed according to the process previously described, using a phase congruency-related function, such as either PC(x) or LE(x).

[0046] Next at 45 a skeletonization step extracts the skeleton of the acquired edge map by removing the non-maximal values of the edge map without breaking the object topology. This is accomplished by setting the non-maximal values to zero in the perpendicular direction (orientation) to the edge. While crystals generally have straight edges, the phase congruency-based edge detection discussed previously finds any kind of edge, not necessarily only straight ones. A line can be distinguished from other objects because the local energy in its perpendicular direction is high, and low in the parallel direction of the line. It is notable that a spot feature has the same local energy for all the orientations. This information is used to remove all the edges that are not line features. Additionally at 46 spur lines may also be removed from the skeleton structure obtained, as a post-edge detection image-cleaning step. In any case, the detected edges are then segmented into discrete, straight line segments in the line segmentation step of 47, followed by a geometric analysis and evaluation at 48, utilizing evaluation criteria characteristic of crystals as previously described. And finally at 49, a determination and decision is made based on the preceding geometric evaluation, as to whether crystals are present in a particular image. It is notable that a second edge detection (not shown) may be optionally utilized in conjunction with the first edge detection process described above, in order to further refine or otherwise limit the edge detection results obtained from the first edge detection.

[0047] While particular operational sequences, materials, temperatures, parameters, and particular embodiments have been described and or illustrated, such are not intended to be limiting. Modifications and changes may become apparent to those skilled in the art, and it is intended that the invention be limited only by the scope of the appended claims.