METHOD FOR REGISTERING SEPARATION PATTERNS
FIELD OF THE INVENTION
The invention relates to methods of registering a pair of one- or two-dimensional separation patterns, for example such patterns as are obtained by one- or two-dimensional electrophoresis.
BACKGROUND OF THE INVENTION
Separation of a solubilized mixture of molecules into its components is routinely perfora ed by forcing the mixture to migrate through a porous solid or fluid medium called the "matrix ". The matrix and other relevant factors are selected so that the various molecules of the -mixture have different migration rates in the matrix according to a particular criterion. The mixture thus separates into its components as it progresses through the matrix. Conditions may be chosen so that the migration rate of a molecule will depend, for example, upon its molecular weight or its electrical charge. At the end of the separation, the components are visualized by various methods either in the matrix or as they are eluted from the matrix, thus producing a one-dimensional separation pattern consisting of an array of spots. Each spot in the array contains all of those molecules in the mixture that migrated through the matrix at a particular rate. The molecules in a given spot may all be identical. More commonly, however, the spot will still contain a mixture of different molecules. If for example the molecules were separated according to molecular weight, then a given spot will contain all the molecules in the original mixture having about the same molecular weight.
Gel electrophoresis is a solid phase separation technique. The matrix, in this case, is a gel slab typically cf polyacrylamide, and the mixture is made to migrate through the gel under the influence of an electric field. In one application, molecules such as proteins or nucleotides separate according to their molecular 5 weight. In another application, lαiown as isoelectric focusing, molecules separate according to their isoelectric point (the pH at which a given molecule is electrically neutral). Other solid phase separation techniques include thin layer chromatography, paper chromatography, and liquid chromatography. Other separation techniques include mass spectroscopy. ι o Since a given spot in a one-dimensional separation pattern may still be a mixture of different molecules, it is common to subject the spots in the pattern to a second round of separation in which the molecules are made to migrate in a direction perpendicular to that of the first round. In the second round, conditions are chosen so that the molecular components of each spot are separated according
1 5 to a different criterion than was used during the first round. After the second round of separation, the molecules are visualized producing a two-dimensional separation pattern. Each spot in the pattern will contain those molecules that migrated through the matrix at the same rates during both separation rounds. In most instances, each spot will contain a single molecular species. 0 For example, a mixture of proteins may first be separated by gel electrophoresis according to their isoelectric point and the resulting one-dimensional separation pattern subjected to a second round of electrophoresis according to molecular weight. More than 5,000 proteins can be resolved on a single gel by this type of two-dimensional electrophoresis. 5 It is often of interest to compare two or more two-dimensional separation patterns. For example, a two-dimensional electrophoretogram of blood plasma proteins from two individuals may be compared. Corresponding spots in the two patterns having different sizes (indicative of a different abundance of a protein in the two samples) or spots having different locations in the patterns (indicative of a
mutated protein) are identified. Such differences may reflect the molecular basis of a disease.
In order to identify altered spots in two separation patterns, the two patterns must first be put in register. This involves identifying pairs of corresponding spots 5 in the two patterns. Generating a one-to-one correspondence between spots in the two patterns (referred to herein as registration) is complicated by the fact that corresponding spots may have different locations and appearances in the two patterns not only due to intrinsic differences in the original samples but also due to unavoidable variability in sample preparation, matrix preparation, running i o conditions, visualization of the pattern and physical warping of the matrix.
Algorithms designed to register a pair of separation patterns are known in the art, for example Melanie™ and PDQuest™ marketed by Bio-Rad Laboratories and Bio Image™ marketed by Bio Image U.K. In these algorithms spots are detected in digital images of the patterns. This yields a "spot-file " for each pattern
15 consisting of a list of the detected spots, their centers, boundaries and volume. The two images are then presented side by side to the user on a monitor screen and the two lists are aligned interactively. The user is prompted to indicate pairs of corresponding spots, which is often a laborious and time-consuming process.
U.S. Patent No. 5,073,963 to Sammons et al. discloses a computerized 0 method for registering spots in a pair of two-dimensional gel electrophoretograms. Their method makes use of coordinate transformation techniques lαiown in the art. A spot file for each pattern is obtained and a transformation of a predetermined type is sought that generates a one-to-one correspondence between the spots in one file and those of the other with minimal distortion of the patterns, as measured by a 5 predetermined distortion criterion.
SUMMARY OF THE INVENTION
Within the context of the present invention, two explicitly described, calculable or measurable variables are considered equivalent to each other when the 0 two variables are proportional to each other.
The present invention provides a method for registering digital images of two separation patterns of the same dimensionality. A pixel is a vector having the same dimensionality as the pattern. The first and second digital images are described by means of gray level function f\ and ?, respectively, defined on a pixel set. that describe the intensity of the image at each pixel. A region of the pixel set is a subset of the pixel set. The appearance of a region in a digital image is the restriction of its gray level function to the region.
In accordance with the invention, the registration is described by means of a transformation T, referred to as the global registration transformation, mapping at least a subset of the pixel set into the pixel set. This is in contrast to prior art methods where registration is performed by applying transformation techniques to the spot lists derived from the digital images. The present invention thus utilizes a larger amount of available data for the registration process than the prior art methods. The utilized data include not only the location and precise structure of the spots of the visualized molecules, but also reproducible staining in the patterns arising from other sources such as streaks and spot shapes.
In accordance with the invention, the global registration transformation T is obtained as follows. At least a portion of the pixel set is divided into m regions Rι,...,Rk,...Rm. For each region Rk, a finite sequence of n transformations SM,...,SA. /,...S/(,»/, is selected, where each transformation S j maps Rk into the pixel set. This generates a sequence of n regions in the pixel set denoted by Skj(Rk) (i=l ,...nk). For each of the n regions S j(Rk) a similarity score is calculated which scores the similarity of the appearance of Rk in the first digital image and the appearance of Sk.j(Rk) in the second digital image. An optimal transformation Sk is selected from among the n transformations St,ι,...,S*, /,...&,«. for which the similarity score of Sk(Rk) is maximal.
The global transformation T is now determined as follows. In each region Rk. a pixel p is selected, pk may be for example the center of gravity of Rk. A global transformation T is then found based upon at least some of the S (pk), for k=l ....m.
The first image may be partitioned into any number m of regions. As extreme examples, it may be partitioned into a single region, or each pixel in the image may be taken as a separate region. The regions may have any shape such as squares, rectangles, circles, ellipses or irregular shapes. The regions in a partition may all be of the same shape, or there may be regions in the partition of different shapes. The regions may or may not overlap. The image may be partitioned into a regular array of identical regions. The partition may be formed taking into account the distribution of spots in the separation pattern. For example, a partition may be formed so that no regions in the partition is devoid of spots while areas in the pattern with many will be divided into a large number of small regions.
In a preferred embodiment, the regions R , are chosen so that each one contains at least two spots, and most preferably three or four spots of the first digital image. If a region Rk contains too many spots, it may not be satisfactorily transformed by any of the transformations S j as assessed by the similarity scores. On the other hand, regions that are have too few spots may have a geometric signature that is not sufficiently unique, in which case several of the transformations S .* may have high similarity scores. For example, a region that contains only one spot can potentially match any other single-spot region. Spots may be detected by methods lαiown in the art. For example, the boundary of a spot in a digital image may be recognized as a location where the Laplacian of the gray level function of the pattern exceeds a threshold value. Thus, for example, the number of spots in the first digital image may be determined and the average area of a region containing four spots calculated. A region of this area is constructed around every spot in the first pattern, and a score calculated for each of the regions that takes into account the number of spots in the region, their contrast and area. For each pair of strongly overlapping regions, the one with the lowest score is discarded.
In accordance with another embodiment, regions in the pixel set are first designated and then ordered. For example, a score is calculated for each of the regions that takes into account the number of spots in the region in the first image,
their contrast and area. A region Ri is selected as the region having a high score near the center of the first digital image. The remaining regions are then sequentially ordered according to increasing distance from the first region.
In yet another embodiment, the transformations Skj are translations of the form
ffk , where « j are translation vectors. The translation vectors may be selected, for example, among lattice points of the form Is, for one-dimensional digital images, and /sι+/zs
2 for two-dimensional images, where s, s* and s
2 are vectors, s* and s
2 are linearly independent, and / an h are integers satisfying |l|<dι and |h|<d for predetermined constants d| and d
2. The vector s or the vectors s* and s
2 are most preferably chosen such that their magnitudes are approximately equal to the radius of the smallest expected spot in the digital image and will typically be around 3 pixels. Other types of transformations S
j are also envisaged within the scope of the invention including linear transformations and affine transformations. The transformations Skj may take into account torn parts or warping of the matrix, for example by using models of deformable membranes as disclosed, for example, in Singh et al., 1998.
The similarity score is designed to take into account differential expression, different staining qualities or techniques, and other irregularities in the two images such as arte factual staining, contamination, air bubbles etc. The scoring function compares the appearance of a region Rk in the first image to that of each of the ii regions Sk.j,(Rι<) in the second image. The scoring functions thus use nearly all of the information contained in the two images. This is in strong contrast with prior art registration methods that use only the list of approximated spot profiles, or Gaussians of the images. The similarity score of Rk in the first image and Skj(Rk) in the second image may be described as a function having, for example, the form ∑G(x, Skj(x)) where the sum is taken over all x e Rk . G is a function that preferably increases if the contrast at x and Sk,*(x), in the first and second images, respectively, both exceed a predetermined minimal contrast. G may be further increased, for example, when the contrast at x and Sk. *(x) are both above a second higher predetermined level. This
favors matching of strong spots over matching of faint spots or spots with dissimilar contrast. The score may also take into account a comparison of the scalar product of the two intensity gradient vectors. This gives an added bonus for spots that have about the same size and shape. This contribution arises from the edge 5 pixels of the spots. In cases where the edges overlap and are similarly oriented, the score is increased. Where there is no correlation or negative correlation, it diminishes. This contribution mimics some of the processing that the human eye does when attempting to match gel images. It is reasonable to prefer a transformation that matches spots with similar contrast, size and shape, to one that i o matches spots that differ widely.
The selection of the optimal transformations Sι,...,Sk,...,Sm may be performed iteratively. For example, once the optimal transformations Sι,...,Sk-ι have been selected, a similarity score may be used that prefers an S that is similar to one or more of any of the Sι,...,Sk-ι.
15 The global transformation T may be obtained from the Sk(pk) by several different methods as are lαiown in the art. If the number of regions (m) is 1, 2 or 3, T may be defined as the unique translation, linear transformation, or affϊne transformation, respectively, determined by the constraints T(pk)=Sk(pk), for k=l,...m. Even when m is large, a translation may be satisfactory, for example, if 0 the two images were obtained from different scans of the same separation pattern. If the images are from different patterns, but the separation technique is lαiown to have very tight geometric tolerances (for example, by using a matrix having a rigid backing), then an affine transformation may be satisfactory. An affine transformation compensates for shift, rotation, scaling, and skew of the axes, but 5 does not compensate for other types of deformation such as warping of the matrix. In cases such as these, T may be a radial basis function transformation, for example, as disclosed in Poggio and Girosi 1989, Poggio and Girosi 1990, or Girosi et al. 1995, or a Delaunay triangulation transformation, as disclosed, for example, in Kanaganathan et al. 1991, or Baker 1987. Once the global
transformation has been found, the image can be transformed in one of several standard lαiown ways such as nearest neighbor, bilinear, or bicubic interpolation.
Once a global transformation function has been obtained in accordance with the invention, it may be revised, for example, using optical flow equations, as
5 disclosed, for example in Aggarwal et al. 1988, or Huang et al, 1994. In this method, a desired correction vector w(x) is found by minimizing the expression
^(Vg • w- dg)1 , wherein Vg* is the gradient vector of the intensity at x and dg is the gray level difference between x and T(x). The minimum is easily found using least mean square fitting techniques for example as disclosed in Press et al., 1992. ι o T(x) is then replaced with T(x)+w(x). This revision improves robustness and precision and may be repeated several times.
The results of the registration may be presented to a user by superimposing the two digital images on a monitor screen such that a pixel in one digital image is superimposed upon its image under the transformation in the other. The two digital
15 images may be false-colored. If, for example, the spots in one image are false colored blue, and those in the other digital image are false-colored yellow, then any superimposed spots in the superimposition will appear green.
The invention may be used to compare each of several separation patterns to one particular separation pattern referred to as the ^master pattern ". The master 0 pattern may be annotated by compiling information on the composition of each spot in the pattern. For example, a two-dimensional electrophoretogram of blood plasma proteins from a healthy individual may be prepared and used as a master pattern. Test electrophoretograms from other individuals are then registered with the master pattern. When the digital image of a test pattern that has been registered 5 with the master pattern is then viewed on a monitor screen, a spot in the test pattern may be selected, and the compiled information on the spot presented to the user.
The invention may be also be used to follow time dependent changes in molecular composition. For example, samples of blood plasma of proteins from an individual undergoing medical treatment may be obtained at various times and a 0 two-dimensional separation pattern of each sample prepared. Each separation
pattern may be registered with a master pattern or the previous pattern in the sequence. The registered patterns may then be displayed sequentially in time on a monitor screen so as to animate time dependent changes in the composition of the samples. It will also be understood that the invention contemplates a computer program being readable by a computer for executing the method of the invention. The invention further contemplates a machine-readable memory tangibly embodying a program of instructions executable by the machine for executing the method of the invention.
Thus, in its first aspect the invention provides a method for registering first and second digital images of first and second separation patterns, respectively, the digital images being described by gray level functions f\ and f∑, respectively, defined on a pixel set, the method comprising the steps of: a dividing at least a portion of the pixel set into m regions Rι,.. .,R
k,...R
m, the m regions each having an appearance in the first image; b For each region Rk, ba selecting a finite sequence of nk transformations
ι,...Sk.m , each transformation mapping R into the pixel set; bb generating nk regions S j(R ) for j=l,...nk in the pixel set, the n
k regions each having an appearance in the second digital image; be calculating a similarity score for each of the n
k regions S
j(R
k) for j=l,...nk, the similarity score scoring the similarity of the appearance of Rk in the first digital image and the appearance
Skj(Rk) in the second digital image; bd selecting a transformation Sk from among the n transformations Su,...,St, ,,...Sk, for which the similarity score of Sk(Rk) is maximal among the similarity scores of the S j(Rk) forj=l,...nk; and
be selecting a pixel pk in Rk; c defining a global transformation T mapping at least a portion of the pixel set into the pixel set based upon at least some of the S (p ); the transformation T having an inverse T" , the second image having an 5 image under T" ; and d obtaining the image of the second image under T" , to produce a modified second image.
In its second aspect, the invention provides a program storage device readable by machine tangibly embodying a program of instructions executable by i o the machine to perform method steps for registering first and second digital images of first and second separation patterns, respectively, the digital images being described by gray level functions f\ and fι, respectively, defined on a pixel set, the registering comprising the steps of: a dividing at least a portion of the pixel set into m regions
15 Ri , ...,Rk,...Rni, the m regions each having an appearance in the first image; b For each region Rk, ba selecting a finite sequence of nk transformations Sk,χ,...,Sk. i,...Sk,M. , each transformation mapping R into the pixel set;
0 bb generating nk regions Skj( k) for j=l,...nk in the pixel set, the n
k regions each having an appearance in the second image; be calculating a similarity score for each of the nk regions Skj(Rk) for j=l,...n
k, the similarity score scoring the similarity of the appearance of R
k in the first image and the appearance 5 Skj(Rk) in the second image; and bd selecting a transformation S from among the n transformations Sk.x,...,Sk, ...Sk,„ι for which the similarity score of S
k(R
k) is maximal among the similarity scores of the
0 be selecting a pixel pk in R
k; c defining a global transformation T mapping at least a portion of the pixel set into the pixel set based upon at least some of the S (p ; the
transformation T having an inverseT-1 , the second digital image having an image under T- 1 ; and d obtaining image of the second image under T
" , to produce a modified second image.
In its third aspect, the invention provides a computer program product comprising a computer useable medium having computer readable program code embodied therein for registering first and second digital images of first and second separation patterns, respectively, the digital images being described by gray level functions f\ and fi, respectively, defined on a pixel set, the registering comprising the steps of: a dividing at least a portion of the pixel set into m regions
Rι,...,Rk,...Rm, the m regions each having an appearance in the first image; b For each region Rk, ba selecting a finite sequence of nk transformations Sk,\,...,Sk, ι,...Sk, , each transformation mapping Rk into the pixel set; bb generating nk regions S j(Rk) for j=l,...nk in the second digital image, nk regions each having an appearance in the second image; be calculating a similarity score for each of the nk regions Skj(Rk) for j=l,...nk, the similarity score scoring the similarity of the appearance of R in the first image and the appearance of
Skj(Rk) in the second image; bd selecting a transformation S
k from among the n transformations Sk,x,...,Sk, ,,...Sk,m for which the similarity score of S (R
k) is maximal among the similarity scores of the
c defining a global transformation T mapping at least a portion of the pixels in the pixel set into the pixel set based upon at least some of the S
k(pk); the transformation T having an inverse T
" , the second digital image having an image under T
" ; and
d obtaining the image of the second image under T
" , to produce a modified second image. In its fourth aspect, the invention provides a program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to perform method steps for registering a sequence of n digital images Iι,...I„ of n separation patterns, respectively, the digital images being described by gray level functions defined on a pixel set, the method comprising the steps of: a registering the digital images l
\ and I
2 according to any one of Claims
1 to 15; to produce an image I'2 that is the modified image of I2; b registering I'k and Ik+i, for all k from 2 to n-1 according to any one of
Claims 1 to 15 to produce an image I'k+i, where I'k and IV+i are the modified images of Ik and Ik+i, respectively. In its fifth aspect, the invention provides a computer program product comprising a computer useable medium having computer readable program code embodied therein for registering a sequence of n digital images Iι,...In of n separation patterns, respectively, the digital images being described by gray level functions defined on a pixel set, the computer program product comprising: a computer readable program code for causing the computer to register the digital images Ii and I2 according to any one of Claims 1 to 15; to produce an image I'2 that is the modified image of I2; b computer readable program code for causing the computer to register
I'k and Ik+i, for all k from 2 to n-1 according to any one of Claims 1 to 15 to produce an image I'k+i, where I'k and I'k+i are the modified images of Iι and Ik+i, respectively. In its sixth aspect, the invention provides a program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to perform method steps for registering a sequence of n digital images
11,... In of n separation patterns, respectively, the digital images being described by gray level functions defined on a pixel set, the method comprising registering a master pattern and the image Ik according to any one of Claims 1 to 15 so as to produce an image I'k, where I'k is the modified image of Ik, for all k from 1 to n.
In its seventh aspect, the invention provides a program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to perform method steps for displaying a sequence of n digital images
I ι....I„ of n separation patterns, respectively, the digital images being described by gray level functions defined on a pixel set, the method comprising: a registering the sequence of images according to Claim 32, so as to produce a second sequence of images I'ι, I'2...I'n, where I'ι is the image Ii and I'k is the modified image of Ik, for k from 2 to n; and b displaying the second sequence of images I'ι,...I'n sequentially on a display.
In its eighth aspect, the invention provides a computer program product comprising a computer useable medium having computer readable program code embodied therein for displaying a sequence of n digital images In,...In of n separation patterns, respectively, the digital images being described by gray level functions defined on a pixel set, the computer program product comprising: a computer readable program code for causing the computer to register the sequence of images according to Claim 32, so as to produce a second sequence of images I'ι, I'2...I'n, where P- is the image Ii and I'k is the modified image of L, for k from 2 to n; and b computer readable program code for causing the computer to display the second sequence of images I'ι,...1',, sequentially on a display. In its ninth aspect, the invention provides, a program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to perfoπΗ method steps for displaying a sequence of n digital images In,... In of n separation patterns, respectively, the digital images being described by gray level functions defined on a pixel set, the method comprising: a registering the sequence of images according to Claim 33 so as to produce a second sequence of images I'ι, I'2,...I'n, where I'k is the modified image of L; and b displaying the second sequence of images I'ι,...I'π sequentially on a display.
In its tenth aspect, t ie invention provides a computer program product comprising a computer useable medium having computer readable program code embodied therein for displaying a sequence of n digital images In,...In of n separation patterns, respectively, the digital images being described by gray level functions defined on a pixel set, the computer program product comprising: a computer readable program code for causing the computer to register the sequence of images according to Claim 33 so as to produce a second sequence of images I'ι, I'2,...I'n, where I'k is the modified image of Ik; and b computer readable program code for causing the computer to display the second sequence of images I'ι,...I'n sequentially on a display.
BRIEF DESCRIPTION OF THE DRAWINGS
In order to understand the invention and to see how it may be carried out in practice, preferred embodiments will now be described, by way of non-limiting examples only, with reference to the accompanying drawings, in which:
Figs, la and lb show a digital image of a first and second two-dimensional separation pattern, respectively, Fig. Id shows a superimposition of the patterns shown in Figs, la and lb, and Fig. lc shows the registration of the two images according to one embodiment of the invention; Figs. 2a and 2b show a digital image of a first and second two-dimensional separation pattern, respectively, Fig. 2d shows a superimposition of the patterns shown in Figs. 2a and 2b, and Fig. 2c shows the registration of the two images according to a second embodiment of the invention; and
Figs. 3a and 3b show a digital image of a first and second two-dimensional separation pattern, respectively, Fig. 3d shows a superimposition of the patterns shown in Figs. 3a and 3b, and Fig. 3c shows the registration of the two images according to a third embodiment of the invention.
DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS
First Embodiment: Global Affine Transformation Model
The first embodiment will be demonstrated in relation to a pair of two-dimensional images. The pixel set is partitioned into a single region (the entire pixel set) and the identity transformation is used as the local transformation S|. The global transformation T is an affine transformation of the form
T is determined as follows. Let x ,...x
π be the n pixels in the pixel set, where x
k -
) . Define To as the identity transformation on the pixel set, and define the gray level function go by go(x )
for every pixel x in the pixel set, where f
\ is the gray level function of the first image. Once T
* and g have been found, T
* +ι and g
*+ι are found as follows. Denote by V(x
L)=(Vι(x),V2(x
k)HV(g
((x))+V(/2(x
k)))/2 the average gradient vector of the two gray level functions g and/2, where fi is the gray level function of the second image. T
*+ι is an affine transformation whose six parameters αu, αι
?2, α
2,ι, aι,ι, bi, and b
2 are found by solving the following overdetermined system of linear equations:
/ . b.χ1 ,(X) xX' χb) xX' Xb) xX' Xx1) V,(x') V-(x i'-) (n J V,(Λ--) χ;V,(x2) Λ-,2V,(λ-:) xy2(x2) V,( 2) V:(r)
Λ-;N,(.v") xX' ,(x") x[Xχx") xX' Xx") V,(λ-") V,(Λ-")
b2
o
*, n is then defined by
A transformation A, is defined as the composition of T
* to T„ A,=T,? T,_ι?- - -?Tι. The process is terminated when an A, is obtained for which the fraction of pixels x in the pixel set for which the distance between x and A,(x ) exceeds a predetermined distance (for example 10 pixels) is small (for example, less than 10% of all of the pixels in the pixel set). The global transformation T is then defined as T=A,.
If at any stage, the transformation T, is large in the sense that regions of the first image are moved large distances in the second image, then the images may be preliminarily blurred and decimated. This may be performed by methods known in the art for example as described in E. Adelson, et al. 1984. The preferred way of blurring and decimating the image is to average every 3X3 subarray of pixels of the images into one pixel of the decimated/blurred image. The image can then be blurred further, by a binomial operator of the form
T, may then be modified by replacing b! and b
2 with 3bι and 3b
2.

Figs, la and 2a schematically show first and second digital images, respectively, of two-dimensional electrophoretograms of urine proteins. The proteins were first separated in the horizontal direction according to their isoelectic point and subsequently in the vertical direction according to their molecular weight. In Figs, la and lb, each spot is represented schematically by a black curve delineating its boundary. The interior of each spot is indicated by hatching that is vertical in Fig. la and horizontal in Fig. lb. Fig. Id shows a superimposition of the images. A global transformation T between the two separation patterns was obtained in accordance with this embodiment and is shown in Fig. lc. In Fig. lc, the appearance of a pixel x in the first digital image has been superimposed on the appearance of its image T(x) in the second image. Regions of overlap of spots in the two images appear in Figs, lc and Id by vertical hatching superimposed on horizontal hatching. The amount of overlap is substantially greater in Fig. lc compared to Fig. Id.
Second Embodiment: Two Dimensional Spline Transformation.
The second embodiment will be demonstrated in relation to a pair of two-dimensional images. The pixel set is partitioned into a grid of rectangles Rk, all having the same dimensions whose centers lie on a lattice spanned by two linearly independent vectors s and s . The local transformations Sk are translations and the global transformation T is a spline.
The width w and height h of the rectangles, can be predetermined constants or determined heuristically in each case. For example, a spot of diameter d (in pixels) in the first image may be selected and the parameters obtained by w = h = \sι\ = \s2\ =10d.
Background contrast images (BCI) Bi and B2 are calculated for the first and second images, respectively. The BCI is a function of the pixels in an image the value of which is 0 for pixels not belonging to a spot. For pixels in a spot, the value of the BCI is a reflection of the contrast of the spot. For two pixels x and y, the pixel score function Psf(x,y) is defined as follows: if Bι(x) or B2(y) is zero, then Psf(x,y) = 0. Otherwise, x is in a spot in the first image and y is in a spot in the second image, and
where (V(/ι(λ-)), V( 2(v))) indicates the absolute value of the scalar product of the two gradient vectors V( ι( )) and V(/χ v)) . This formula allows a fixed contribution (1) to the score if x and y both belong to a spot. This allows spots representing differently expressed proteins to match. There is also a "bonus" for cases where both spots have high contrast, by increasing the score in correlation with the minimum of the two contrast values. This encourages matching of strong spots over matching of faint spots or spots with dissimilar contrast. There is also a contribution from the scalar product of the two gradient vectors. This gives an added bonus for spots that have the same size and shape.
Where there is no correlation or negative correlation, the score is low. This mimics some of the processing that the human eye does when attempting to match images. The constants a, β are weighting factors that determine the relative contribution of the three summands in the definition of Psf. a will typically be set to the value of the minimal contrast required for a spot. Typical values of α and β are a = β = 0.1 .
The local transformations Skj are translations of the form
for j= l ,...nk. The « ,j are selected so that the centers of the regions Skj(Rk) in the second image form a lattice with a typical distance of s pixels between neighbors, where s is the radius of the smallest expected spot on the gel. This decreases the computational load by s
2 (typically, s=3). The «kj are also selected so that the regions S
kj(Rk) cover an area large enough that S (Rk) lies far from the boundary of the union of all of the regions Skj(Rk). The score function for S j is the sum of the S
kj(x)) for all x in R
k. An S j having maximal score is then selected and designated as S . As in the first embodiment, this selection of Sk may be refined, for example, by optical flow techniques as disclosed, for example, in Aggarwal et al., 1988, Fluang et al. 1994. Alternatively, the S may be refined taking into account the scores of rectangles in the neighborhood of the selected one and determining optimal shift. Once an S
k is obtained, it may be checked. For example, if it is on the boundary of the original search area, or its score is below a predetermined threshold, or it deviates substantially from its nearest neighbors, it may be discarded.
The global transformation in this embodiment is a two dimensional spline calculated from the transformations Sk by methods lαiown in the art, for example as disclosed in de Boor, 1978. Since the rt are computed along a grid, the spline can be computed first on the horizontal lines, and than on the vertical lines.
Figs. 2a and 2b schematically show first and second digital images, respectively, of two-dimensional electrophoretograms of urine proteins. The proteins were first separated in the horizontal direction according to their
isoelectic point and subsequently in the vertical direction according to their molecular weight. In Figs. 2a and 2b, each spot is represented schematically by a black curve delineating its boundary. The interior of each spot is indicated by hatching that is vertical in Fig. 2a and horizontal in Fig. 2b. Fig 2d shows a superimposition of the images. The first separation pattern has been divided into rectangles whose centers lie on a lattice. A global transformation T between the two images was obtained in accordance with this embodiment and is shown in Fig. 2c. In Fig. 2c, a pixel x in the first digital image has been superimposed on its image T(x) in the second image. Fig. 2c also shows the system of rectangles that was used in this example. Regions of overlap of spots in the two images appear in Figs. 2c and 2d by vertical hatching superimposed on horizontal hatching. The amount of overlap is substantially greater in Fig. 2c compared to Fig. 2d.
Third Embodiment: Delaunay Triangulation Transformation Model
The third embodiment will now be described in relation to a pair of two-dimensional images. As in the second embodiment, in this embodiment the pixel set is partitioned into rectangles R . However, the rectangles are selected in accordance with the data and do not necessarily form a grid. The local transformations Sk., are translations and the global transformation T is a spline.
The rectangle sequence determination proceeds as follows: A spot list is first prepared of the first image by methods lαiown in the art, for example as described in Gonzalez et al., 1977 or Pratt, 1991. The average area S of a rectangle (in pixels) containing four spots is S = 4A / n , where A is the number of pixels in the pixel set, and n is the total number of spots in the image. A rectangle is constructed with area A around every spot. In each rectangle, all spots are found with contrast above a given threshold. A score s is calculated for the rectangle by s = «[/?] • c , where n is the number of spots in the rectangle, c is the average contrast of the spots, and a[ιx] is a weighting function that prefers rectangles containing a number of spots close to predetermined number. In a
preferred embodiment, rectangles are preferred in which the number of spots is around 3. In this case, a can be defined, for example, by 0, 10, 15,20, 15, 12, 10, 8, 7. 6, 5, 4, 3, 2, 1, for n equal to 0 to 14, and α[n]=0 for n greater than 14. For each pair of strongly overlapping rectangles, the one with the higher score is retained and the one with the lower score is deleted.
The rectangles are now ordered. A rectangle Ri with a high score is found near the center of the image. The remaining rectangles are subsequently ordered by increasing distance from Rj.
The global transformation in this embodiment is an affine transformation and a Delaunay interpolation triangulation transformation. The affine transformation is obtained as in the first embodiment, and then the Delaunay interpolation is used for the differences between the Sk and the affine transformation. This provides a stable and robust transformation that is easy extrapolated, as well as a method for dealing with relatively small and non-uniform perturbations. The global transformation in each of the Rk is obtained in a process similar to that of the first embodiment. However, in this embodiment, the global transformation is recomputed after each Sk has been found. Since initially there may be insufficient parameters to compute the affine transformation, simpler transformations are computed in the first three steps (identity, translation, linear transformation). At each stage, Sk is used in obtaining the shift vector «k+ι of Sk+i
A Delaunay triangulation is obtained by methods known in the art, for example, as described in Kanaganathan et al, 1991. The Delaunay triangulation is used to interpolate the global transformation T in the convex hull of the triangulation vertices. At the four corners of the image, the Delaunay triangulation is set to zero so that the corners of the image, the global transformation T is affine only.
Figs. 3a and 3b schematically show first and second digital images, respectively, of two-dimensional electrophoretograms of urine proteins. The proteins were first separated in the horizontal direction according to their isoelectic point and subsequently in the vertical direction according to their
molecular weight. In Figs. 3a and 3b, each spot is represented schematically by a black curve delineating its boundary. Fig. 3d shows a superimposition of the images. The interior of each spot is indicated by hatching that is vertical in Fig. 3a and horizontal in Fig. 3b. The pixel set was divided into rectangles of different sizes so that each rectangle has about the same number of spots. A global transformation T was obtained in accordance with this embodiment and is shown in Fig. 3c. In Fig. 3c, the appearance of a pixel x in the first digital image has been superimposed on the appearance of its image T(x) in the second image. Fig. 3 also shows the system of rectangles that was used in this example. Regions of overlap of spots in the two images appear in Figs. 3c and 3d by vertical hatching superimposed on horizontal hatching. The amount of overlaps is substantially greater in Fig. 3c compared to 3d.
In the method claims that follow, alphabetic characters used to designate claim steps are provided for convenience only and do not imply any particular order of performing the steps.
References
Adelson, E. et al. RCA Engineer 29, No.6, Nov/Dec 1984, pp 33-41.
Aggarwal .1. K., and Nandhakumar N., "On the computation of motion from sequences of images - A review", Proceedings of the IEEE, 76(8), pp. 917-935, 1988.
Baker T.J., "Three dimensional mesh generation by triangulation of arbitrary point sets" AIAA 8th Computational Fluid Dynamics Conference 1987.
de Boor. C, "A Practical Guide to Splines", Springer- Verlag, New- York, NY, 1978.
Girosi, F., Jones, M. and Poggio, T., "Regularization Theory and Neural Networks Architectures", Neural Computation, Vol. 7, 219-269, 1995. Gonzalez. R.C.and Woods, R.E., Digital Image Processing, Addison-Wesley, 1992.
Gonzalez. R.C. and Wintz, P.A.. "Digital Image Processing", Reading, MA., Addison- Wesley 1977.
Huang, T.S.. and Netravali, A.N., "Motion and Structure from Feature Correspondences: A Review", Proceedings of the IEEE, 82(2), pp. 252-268, 1994.
Kanaganathan S., and Goldstein N.B., "Comparison of four point adding algorithms for Delaunay type three dimensional mesh generators" , IEEE Transactions on magnetics, Vol 27, No 3, May 1991.
Poggio T. and Girosi, F, "A Theory of Networks for Approximation and Learning", Al memo 1 140, July 1989.
Poggio, T. and Girosi, F., "Networks for Approximation and Learning", Proceedings of the IEEE, vol. 78, no. 9, pp. 1481-1497, September 1990.
Press W.H.. Flannery B.P., Teukolsky S.A., and Vetterling W.T., "Numerical Recipes In C (Second Edition) ", Cambridge University Press, Cambridge, 1992.
Pratt, W.K., "Digital Image Processing (Second Edition) ", New York, Wiley 1991.
Singh, A., Goldgof, D., and Teraopoulos, D., "Deformable Models in Medical Image Analysis", IEEE Computer Society, 1998.