IMAGE ANALYSIS SYSTEM AND METHOD

Abstract

An image analysis system includes a video camera that collects YUV color images of a liquid sample disposed between a capital and a pedestal, the color images being collected while a light source shines light through an optical beam path between the capital and the pedestal, and a processor adapted to i) obtain from the YUV color images a grayscale component image and a light scatter component image, and ii) obtain at least one binary image of the grayscale component image and at least one binary image of the light scatter component image.

Claims

1.-24. (canceled)

25. A spectrometer including an image analysis system, the system comprising: a. a video camera that collects images of a liquid sample disposed between a capital and a pedestal of the spectrometer, the images being collected while a light source shines light through an optical beam path between the capital and the pedestal for photometric or spectrometric measurement; and b. a processor adapted to detect any bubble in a column of the liquid sample using the images.

26. The image analysis system of claim 25, wherein the images are YUV color images, and the processor is further adapted to i) obtain from the YUV color images a grayscale component image and a light scatter component image, and ii) obtain at least one binary image of the grayscale component image and at least one binary image of the light scatter component image.

27. The image analysis system of claim 26, wherein the at least one binary image of the grayscale component image includes a first binary image of the grayscale component image obtained from applying an upper dynamic threshold and a lower dynamic threshold obtained from an interpolation between left and right background variation thresholds in the grayscale component image.

28. The image analysis system of claim 27, wherein the at least one binary image of the grayscale component image includes a second binary image of a grayscale isotropic gradient image, the grayscale isotropic gradient image obtained using the grayscale component image and a static threshold based on isotropic gradient image background noise statistics.

29. The image analysis system of claim 28, wherein the at least one binary image of the grayscale component image includes a composite binary image obtained from a combination of the first and second binary images.

30. The image analysis system of claim 29, wherein the processor is further adapted to detect location of the column of the liquid sample and location of the optical beam path from the composite binary image.

31. The image analysis system of claim 30, wherein the processor is further adapted to detect any bubble in the column of the liquid sample using both the grayscale component image and the at least one binary image of the light scatter component image.

32. The image analysis system of claim 31, wherein using the grayscale component image includes applying a ring detection filter to a grayscale isotropic gradient image obtained from the grayscale component image.

33. The image analysis system of claim 31, wherein using the at least one binary light scatter component image includes applying a morphological filter to the at least one binary image of the light scatter component image.

34. The image analysis system of claim 31, wherein the processor is further adapted to distinguish a bubble in the optical beam path from a bubble out of the optical beam path using the grayscale component image, the at least one binary image of the light scatter component image, and the calculated location of the optical beam path.

35. A method of analyzing an image, the method comprising: a. collecting images of a liquid sample disposed between a capital and a pedestal of a spectrometer, the images being collected while a light source shines light through an optical beam path between the capital and the pedestal for photometric or spectrometric measurement; b. detecting location of a column of the liquid sample and location of the optical beam path from the images; c. detecting any bubble in the column of the liquid sample using the images; and d. reporting an image analysis summary to a display.

36. The method of analyzing an image of claim 35, wherein collecting images of the liquid sample includes collecting YUV color images, obtaining from the YUV color images a grayscale component image and a light scatter component image, and obtaining at least one binary image of the grayscale component image and at least one binary image of the light scatter component image.

37. The method of analyzing an image of claim 36, wherein the at least one binary image of the grayscale component image includes a first binary image of the grayscale component image obtained from applying an upper dynamic threshold and a lower dynamic threshold obtained from an interpolation between left and right background variation thresholds in the grayscale component image.

38. The method of analyzing an image of claim 37, wherein the at least one binary image of the grayscale component image includes a second binary image of a grayscale isotropic gradient image, the grayscale isotropic gradient image obtained using the grayscale component image and a static threshold based on isotropic gradient image background noise statistics.

39. The method of analyzing an image of claim 38, wherein the at least one binary image of the grayscale component image includes a composite binary image obtained from a combination of the first and second binary images.

40. The method of analyzing an image of claim 39, further including detecting location of a column of the liquid sample and location of the optical beam path from the composite binary image.

41. The method of analyzing an image of claim 40, wherein detecting any bubble in the column of the liquid sample includes using both the grayscale component image and the at least one binary image of the light scatter component image.

42. The method of analyzing an image of claim 41, wherein using the grayscale component image includes applying a ring detection filter to a grayscale isotropic gradient image obtained from the grayscale component image.

43. The method of analyzing an image of claim 41, wherein using the at least one binary light scatter component image includes applying a morphological filter to the at least one binary image of the light scatter component image.

44. The method of analyzing an image of claim 41, further including distinguishing a bubble in the optical beam path from a bubble out of the optical beam path using the grayscale component image, the at least one binary image of the light scatter component image, and the calculated location of the optical beam path.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0015] FIG. 1A is an illustration of a cutaway section of the optical path in a prior art spectrophotometer.

[0016] FIG. 1B is a schematic illustration of a spectrophotometer where a sample in a form of a liquid drop (column) is contained in a compartment between two planar interface surfaces. The source of a UV light is located on the top part (capital) and a spectrophotometric detector is beneath the bottom part (pedestal) of the device. An optical path can be established between the source and the detector. A video camera collects color images of the liquid drop against a background.

[0017] FIG. 2 shows original high-resolution camera image of the measurement compartment. The indicated rectangular area shows the initial region-of-interest (ROI), which is large enough to capture all the necessary features of the image, including a large portion of its gray background.

[0018] FIGS. 3A-3C show examples of normal liquid columns that provide good optical paths between the two interface surfaces (capital and pedestal).

[0019] FIGS. 4A-4C show examples of a first type of possible defects of the liquid columnskewed (FIG. 4A) and off-center (FIG. 4C) columns that provide only partial or no optical path between the two interface surfaces.

[0020] FIGS. 5A-5C show examples of a second type of defectsscattered light that is caused by the reflection from the surface of the bubbles or other inclusions. Bubbles are made of air (or gas) and are considered defects that degrade spectrophotometric measurement.

[0021] FIGS. 6A-6C show examples of two terminal conditionsa broken column (split apart in two drops) (FIGS. 6A and 6B) or void/empty (FIG. 6C) space between the two interface surfaces. No spectrophotometric measurement can be made in these conditions.

[0022] FIGS. 7A-7C illustrate the component image extraction that is used for the image analysis. FIG. 7A is the ROI part of the camera source color image. FIG. 7B is the grayscale (intensity) component that is extracted from FIG. 7A. FIG. 7C is its extracted complementary light scatter component (used for the scattered light measurement).

[0023] FIG. 8A is a flowchart of image collection from the camera video stream and grayscale and light scatter component image extraction from the YUV color images.

[0024] FIG. 8B is a flowchart of the column image analysis method.

[0025] FIG. 9 shows the extracted initial ROI image that contains the grayscale (luma) component of the original camera image.

[0026] FIG. 10 shows the histogram of the initial grayscale ROI image (from FIG. 9). The largest peak represents background intensity. The delineated area between T.sub.min and T.sub.max shows the selected background intensity range.

[0027] FIG. 11 illustrates the background area on the grayscale image. Background pixels are marked by stippling using the thresholds that match the selected range. Note that background intensity range is underestimated to prevent the features-of-interest from being marked as background.

[0028] FIG. 12 illustrates how the horizontal part of the ROI with the features gets extracted using the cumulative profile from the grayscale image (FIG. 11). The profile is made by summation of the pixel values for each vertical column of the image. The background value is being subtracted from the pixel values in order to obtain a good ratio between valleys and peaks on the profile plot. The target area represents the raised central plateau in the profile.

[0029] FIG. 13 shows the extracted horizontal part of the grayscale image that is used for processing and further extracting the ROI of the compartment features such as the capital and the pedestal.

[0030] FIG. 14 represents the result of applying horizontal gradient filter to the grayscale image (from FIG. 13). It highlights two main edgesone is for the capital (top horizontal edge 1420) and second for the pedestal (bottom horizontal edge 1410).

[0031] FIG. 15 shows how the vertical part of the ROI gets extracted using the profile of the cumulative profile from the grayscale image (FIG. 11). The profile is made by summation of the pixel values for each horizontal column of the image. The vertical portion of the target ROI resides between the two tallest peaks, which represent the horizontal edges of the capital and the pedestal.

[0032] FIG. 16 shows the extracted grayscale ROI image that is used for further extraction of the instrument interface area.

[0033] FIG. 17 shows the resulting image after applying the diagonal (45 degree) gradient filter and shows the extracted right diagonal edge of the upper conical capital surface that is fitted by a line segment (shown in white).

[0034] FIG. 18 shows the resulting image after applying the diagonal (135 degree) gradient filter and shows the extracted right diagonal edge of the left conical capital surface that is fitted by a line segment (shown in white). The distance (in pixels) between the two ends of the extracted edge segments serves as a measure of the capital diameter 1110 shown in FIG. 11. The measured capital diameter 1110 is then used for accurate extraction of the sample compartment area.

[0035] FIG. 19 shows the final extracted ROI grayscale image that is centered (relatively to the axis of the light path) and contains only the area where the liquid column (droplet) is expected with a sufficient amount of the background area on both sides of the image.

[0036] FIG. 20 illustrates the binarization (by a thresholding) of a column grayscale image. It shows the rectangular areas (white rectangles) on both horizontal sides of the image where the background parameters are being calculated. Those parameters are used for setting the binarization thresholds. Note that the result highlighted stippled area (that represents the binary image foreground) has a lot of disconnected segments. It is hard to connect those segments into a column shape.

[0037] FIG. 21 illustrates the dynamic thresholding technique for just one row of grayscale image. The threshold band is created by using mean values and the cumulative standard deviation () of pixel intensity values from two background regions (left and right, marked in gray) and applying a linear interpolation between upper and lower band thresholds. The band width is 6. Pixels with intensity above (and below) the band are assigned to a binary foreground. Note the difference in image brightness between the left and right parts of the image.

[0038] FIG. 22 shows a gradient image that is obtained by applying an isotropic gradient filter. The standard deviation of the gradient intensities is calculated using the background areas, which is used for calculating a threshold for the image binarization. Every pixel whose intensity is above the calculated threshold is assigned to a binary foreground (which is shown as a stippled overlay on top of the gradient image).

[0039] FIG. 23 shows the result of combining two binary foreground images, the one from the grayscale image and another one from the gradient image (FIG. 22) to produce a composite binary image. The resulting binary image is improvedit is much easier now to restore the column shape.

[0040] FIG. 24 shows the composite binary image overlaid on top of the grayscale image with the additional application of two horizontal foreground stripes on the top and bottom borders of the binary image. It enables further connecting the left and right edges of the column.

[0041] FIG. 25 shows the result of applying the hole filling operator to the binary image and removing the artificial stripes that are close to the borders. Subsequent morphological filtration removes residuals and smooths the foreground edges.

[0042] FIG. 26 shows the result of fitting a rectangle into the column shape that represents the theoretical optical path (a cross-section of a cylinder in 3D space).

[0043] FIG. 27 shows an example of a defect in a liquid column with a distinctive bubble on the top of the image.

[0044] FIG. 28 shows a gradient image of the grayscale image with a bubble on the top that becomes a ring-shaped feature.

[0045] FIGS. 29A-29B illustrate the ring detection filter kernel (size of 99 pixels). FIG. 29A illustrates a 3D geometric shape of the filter while FIG. 29B is its digital approximation as a 99 matrix (for illustrative purposes the matrix values are multiplied by 100). The 3D shape of the filter represents a raised circular plateau around its perimeter and a negative circular cavity area in the center. The rest of the kernel is set to the intensity value of 0. The filter is designed to correlate with ring-like particles on the intensity gradient image (see FIG. 27) that enables measuring the ring-shaped feature intensity.

[0046] FIGS. 30A-1 through 30B-2 illustrate the application of the ring filter to a single image row: the gradient image (FIG. 30A-1) with its row profile (FIG. 30A-2) is shown in the top left. The ring filter (top right) is applied to all row pixels by using the formula shown below, where K.sub.j is a kernel coefficient in kernel position j, P.sub.i is a pixel intensity in row position i and F.sub.i is a result value in position i. The ring filter is applied to all image rows. FIG. 30B-1 shows the resulting image with its cross-section row. Note the result peak, shown in FIG. 30B-2, that corresponds to the ring-shape feature on the gradient image (marked with an arrow on top).

[0047] FIG. 31 illustrates obtaining a light scatter binary image from the blue chroma component image. First, the binarization threshold values are calculated using the background area standard deviation and the image gets binarized using the threshold. The resulting binary image is cleared by using the morphological closing filter followed by removing small foreground features.

[0048] FIGS. 32A-32B illustrate the case of a broken column. FIG. 32A is an initial extracted ROI image. FIG. 32B is a (binary) image that shows the result of the analysis algorithm application that correctly detects two broken pieces (highlighted by stippling).

[0049] FIG. 33 is an exemplary image of a liquid column.

[0050] FIG. 34 is a resulting gradient image of the image shown in FIG. 33.

[0051] FIG. 35 is a resulting binary image of the image shown in FIG. 33.

[0052] FIG. 36 is an exemplary image of a liquid column with bubbles.

[0053] Like reference numerals refer to corresponding parts throughout the several views of the drawings.

DETAILED DESCRIPTION OF EMBODIMENTS

[0054] In the description of the invention herein, it is understood that a word appearing in the singular encompasses its plural counterpart, and a word appearing in the plural encompasses its singular counterpart, unless implicitly or explicitly understood or stated otherwise. Furthermore, it is understood that for any given component or embodiment described herein, any of the possible candidates or alternatives listed for that component may generally be used individually or in combination with one another, unless implicitly or explicitly understood or stated otherwise. Moreover, it is to be appreciated that the figures, as shown herein, are not necessarily drawn to scale, wherein some of the elements may be drawn merely for clarity of the invention. Also, reference numerals may be repeated among the various figures to show corresponding or analogous elements. Additionally, it will be understood that any list of such candidates or alternatives is merely illustrative, not limiting, unless implicitly or explicitly understood or stated otherwise. In addition, unless otherwise indicated, numbers expressing quantities of ingredients, constituents, reaction conditions and so forth used in the specification and claims are to be understood as being modified by the term about.

[0055] Accordingly, unless indicated to the contrary, the numerical parameters set forth in the specification and attached claims are approximations that may vary depending upon the desired properties sought to be obtained by the subject matter presented herein. At the very least, and not as an attempt to limit the application of the doctrine of equivalents to the scope of the claims, each numerical parameter should at least be construed in light of the number of reported significant digits and by applying ordinary rounding techniques. Notwithstanding that the numerical ranges and parameters setting forth the broad scope of the subject matter presented herein are approximations, the numerical values set forth in the specific examples are reported as precisely as possible. Any numerical values, however, inherently contain certain errors necessarily resulting from the standard deviation found in their respective testing measurements.

[0056] As the popularity of UV/Vis spectrometers such as the NanoDrop (Thermo Electron Scientific Instruments, Madison Wis.) grows, there is a demand to improve the reliability of its measurement technique that relies, as discussed above, on the surface tension of a liquid sample (drop). A problem is that the spectral measurement of the liquid drop (column) can be degraded by unevenness of the liquid column shape and its position during the measurement. The liquid column can be misshaped (skewed), off-center (relative to the light path of the instrument), can contain bubbles or other inclusions, or the liquid column can even be broken apart. Presently available instruments have no ability to automatically identify and report these problems, while the visual inspection of the liquid drop shape is very limited and unreliable. The design described herein includes a high-resolution video camera in front of the sample compartment that focuses on the liquid column and uses computer vision algorithms for automatic identification of the liquid column defects and reporting them to the instrument operator. The image analysis system also includes measurement and reporting of scattered light that is caused by bubbled liquid, which degrades the measurement. The video camera and spectrometer are synchronized together and every spectral measurement is accompanied with a column image quality metric.

[0057] FIG. 1B shows a spectrophotometer 100 where a sample in a form of a liquid drop (column 9) is contained between two planar interface surfaces 2 and 7 (compartment). The source 10 of light is located on the top part (capital 2) and a spectrophotometric detector 20 is beneath the bottom part (pedestal 4) of the device and an optical path 30 can be established between them.

[0058] The quality of the measurement depends on the quality of the measured liquid column 9 during the time when the light beam 3 is passing through it. The quality of the column 9 is very hard to analyze visually because the gap (distance between two interface surfaces 2 and 7) is too narrow1 mm or less.

[0059] The possible column defects can be summarized in the following categories:

[0060] Skewed and off-centered column. Examples of this defect are presented in FIGS. 4A-4C;

[0061] Bubbled column of column with inclusions. See FIGS. 5A-5C for examples.

[0062] Broken column or empty compartment (no liquid drop). This is a terminal defectno column quality measurement will be made. See FIGS. 6A-6C for examples.

[0063] Turning back to FIG. 1B, a video camera 40 collects YUV color images of a liquid sample 9 disposed between a capital 2 and a pedestal 4, the color images being collected while a light source 10 shines light 3 across an optical beam path 30 between the capital 2 and the pedestal 4 for inspection of the quality of the liquid column 9. The camera 40 is mounted on the instrument pedestal 4 and focused into the gap between the two light transmitting interfaces 2 and 7 where the liquid column 9 is expected to be formed. FIG. 2 shows a typical video camera image. Turning back to FIG. 1B, the camera view is shielded from the direct ambient light by a cover 50 that also creates a background 60. The surface of the background 60 is made of a homogeneous reflective material and has a light (about 82% brightness and 12% gloss) gray color. The background surface 60 is similar to the gray card that is used in photography for creating diffuse (Lambertian) reflectance. The background surface light gray homogeneous reflective background minimizes glossiness and makes the background homogeneous, which produces a narrow, Gaussian-shaped peak on the intensity histogram when the image is acquired by a camera.

[0064] The camera video can be inspected by the operator of the instrument; however, a more accurate and convenient way is to use machine vision algorithms. A processor 70 is adapted to i) obtain from the YUV color images a grayscale component image and a light scatter component image, and ii) obtain at least one binary image of the grayscale component image and at least one binary image of the light scatter component image.

[0065] Scattered light is emitted while acquiring a spectrum on a column with inclusions, such as bubbles of gas or air. The blue wavelength range is the prevailing component in the scattered light, because of the strong (.sup.4) wavelength dependence of Rayleigh scattered light (shorter wavelengths (blue) are Rayleigh scattered more strongly than longer wavelengths (red)). The resulting spectrum quality can be degraded due to the loss of the beam energy that reflects from the bubbles and gets scattered. By applying machine vision algorithms, it is possible to quantitatively measure the amount of scattered light.

[0066] Although it is possible to analyze the original color image on a modern computer, it leads to unnecessary complexity and redundancy. For image analysis, two intensity-only (grayscale) images are created in one of three possible ways:

[0067] 1. When one has just one color RGB image (snapshot), the image is extracted and the following two component images are created:

[0068] a. A grayscale component (luma) image (L) is created by averaging the Red (R) and Green (G) components from the original RGB image. For every x,y-positioned image pixel, the following calculation is applied:

L(x,y)=(R(x,y)+G(x,y))/2;

[0069] b. The blue chromatic component (for the light scatter component) image (S) is created by using the original blue (B) component from the RGB image and calculating the following complementary image, for every x,y-positioned pixel as follows:

S(x,y)=max(0,B(x,y)L(x,y));

[0070] 2. In the case of the YUV image format (that is available on the Android/Linux system) the calculation of the two component images is:

[0071] a. The grayscale component image is the Y (luma) component of the original YUV image, that is: L(x,y)=Y(x,y);

[0072] b. The light scatter component image S is created by using the U-chrominance component from the YUV image and calculating the following complementary image, for every x,y-positioned pixel: S(x,y)=max(0, U(x,y)128);

[0073] 3. In case a sequence of YUV images of varying flash light is available from the camera video stream obtained according to the flowchart shown in FIG. 8A (step 800 in FIG. 8B), 3 images for instance, the calculation of the two component images is done as follows:

[0074] a. A grayscale component image is calculated as an average of all available Y.sub.i (luma) components of the original YUV images from the sequence (step 801 in FIG. 8B), that is:

L(x,y)=(Y.sub.1(x,y)+Y.sub.2(x,y)+Y.sub.3(x,y))/3;

[0075] b. Let S.sub.i(x,y) be the light scatter (blue chromatic) intensity for the i-image in pixel (x,y) that is calculated from the U.sub.i-component using the formula above (see 2.b). Then the S(x,y) light scatter component image is calculated by taking the maximum of all available S.sub.i(x,y) for each pixel (x,y) as follows:

S(x,y)=max(S.sub.1(x,y),S.sub.2(x,y),S.sub.3(x,y));

[0076] The maximum of S.sub.i(x,y) is used to obtain the maximum scattered light that corresponds to the flash occurrence moment.

[0077] See FIGS. 7B and 7C for illustrations of the two-component split images. FIG. 7A is the fragment of the source camera color image. FIG. 7B is the grayscale component image that is extracted from the source RGB image and FIG. 7C is its light scatter component.

[0078] The following steps form the liquid column analysis algorithm:

[0079] 1. Let a grayscale image L of size MN consist of pixels g.sub.i,j, such that i[0, M1], j[0, N1]) 0g.sub.i,j255.

[0080] In other words g.sub.i,j is the image rectangular area and its value (intensity) can vary from 0 to 255.

[0081] 2. Use the extracted grayscale image L (FIG. 9) to calculate a histogram shown in FIG. 10 of all pixels of the source image. Find the image background intensity range (T.sub.min, T.sub.max) using the largest intensity peak of the image histogram (FIG. 10). The background is assumed to be homogeneous with a Gaussian intensity distribution and the background area is assumed to take at least 50% of the extracted grey image. Check the histogram parameters: peak position and its width to verify the background quality. If the peak position is less than 40 (too dark) or greater than 210 (too bright) or the peak width is wider than 80, then report that the background is poor and skip the rest of image analysis. (steps 803 and 804 in FIG. 8B).

[0082] 3. Create a horizontal cumulative profile by summation of the absolute pixel-minus-background values along each image column. That is, by calculating .sub.i.sub.i,j where

[00001]${}_{i,j}=\{\begin{array}{c}{g}_{i,j}-T_{\mathrm{ma}.Math..Math.x},\mathrm{if}.Math..Math.g_{i,j}>T_{\mathrm{ma}.Math..Math.x}\\ T_{m.Math..Math.i.Math..Math.n}-g_{i,j},\mathrm{if}.Math..Math.g_{i,j}<T_{m.Math..Math.i.Math..Math.n}\\ 0,\mathrm{otherwise}\end{array};$

[0083] 4. See FIG. 12 for an example of a cumulative profile that is calculated for the image shown in FIG. 11.

[0084] 5. Find the left and right edges of the instrument capital (upper part, see FIG. 11) by finding edges of the plateau intensity in the profile (see FIG. 12). They correspond to the capital position and its diameter on the grayscale image and are used for the region-of-interest (ROI) extraction.

[0085] 6. Extract a new ROI image that embraces the found features (the capital and the pedestal) with additional extension on both sides as shown in FIG. 13.

[0086] 7. Apply a horizontal gradient (Sobel operator, see page 578 of Digital Image Processing, Rafael C. Gonzalez and Richard E. Woods, 2.sup.nd Ed., Prentice Hall, 2002, (hereinafter Gonzalez) the entire contents and teachings of which are hereby incorporated by reference in their entirety) filter to the extracted ROI grayscale image to find the lower horizontal capital edge 1420 and the upper horizontal pedestal edge 1410 thereby detecting a region-of-interest that includes a location of the capital and the pedestal from the grayscale component image as shown in FIG. 14.

[0087] 8. Create a vertical cumulative profile shown in FIG. 15 of the gradient image shown in FIG. 14 by summation of pixel values along each image row, that is .sub.ig.sub.i,j, where g.sub.i,j is a pixel intensity of the gradient image at the i,j position.

[0088] 9. Find two main intensity peaks 1510 and 1520 on the vertical cumulative profile (FIG. 15). Its right (tallest) peak 1520 corresponds the lower edge of the capital (top boundary or top horizontal edge) while the left (second intensity) peak 1510 corresponds to the upper edge of the pedestal (bottom boundary or bottom horizontal edge).

[0089] 10. Using the found boundaries from the previous step, extract the vertical part of the ROI sub-image for further processing (FIG. 16).

[0090] 11. Find the right edge of the capital on the top of the ROI image (FIG. 16). Apply a 45-degree diagonal gradient filter on the image. Fit the top right set of the gradient intensity pixels with a diagonal line segment by finding the best least-square fit to the right diagonal edge 1710 (FIG. 17).

[0091] 12. Continue on to finding the left edge of the capital. Apply a 135-degree diagonal gradient filter. Fit the top left set of the gradient intensity pixels with a 135-degree diagonal line segment by finding the best least-square fit. FIG. 18 shows both the fitted capital left diagonal edge 1810 and the right diagonal edge 1820. Using the left and right ends of the found segments, calculate the capital diameter in pixels 1830 and X.sub.C the position of the optical path center 1840. The calculated capital diameter 1830 enables accurate matching of the real instrument coordinate units with the image (pixel) units. The optical path center X.sub.C is used for the column feature position evaluation.

[0092] 13. Use previously found capital diagonal edges and knowledge of the actual sizes of the capital and pedestal to extract the final ROI image that is centered with respect to both the instrument capital and pedestal (step 805 in FIG. 8B). Add extra room on the left and the right of the extracted ROI image in order to include more background area as shown in FIG. 19.

[0093] 14. Calculate background parameters for the left and the right parts of the image using the left background rectangular areas 2010 and 2020 and right background rectangular areas 2030 and 2040 where the background is expected (see white rectangles in FIG. 20). The selected background rectangular areas enable estimating the background intensity ranges for the different parts of the grayscale image, while accounting for illumination differences between left and right parts of the image.

[0094] 15. A thresholding technique is then applied to create convenient binary images (an example is shown as a stippled image overlay in FIG. 20), which are used for the column positioning and measurement. A binary image is a digital image that has only two possible values for each pixel. Usually those two values are called foreground (value=1) and background (value=0). See Gonzalez, page 595, for a description of thresholding techniques.

[0095] 16. Create a first binary image of the grayscale ROI image by applying dynamic thresholding that uses an interpolation between the left and right background thresholds to calculate a threshold for each pixel individually. FIG. 21 illustrates the dynamic thresholding technique for just one row 2110 of grayscale image. The threshold band is created by using mean values and the overall standard deviation () of pixel intensity values from two background regions 2120 and 2130 (left and right rectangles, respectively, that are marked in gray) and applying an interpolation between the upper and lower band thresholds. The band width is 6 . Pixels with intensity above (and below) the band are assigned to a binary foreground. Note the difference in image brightness between left and right parts of the image of FIG. 21. The left lower threshold value T.sub.min. 2150 is calculated by using the left background area mean value M.sub.L 2145 and its standard deviation as T.sub.min=M.sub.L3. The left upper threshold T.sub.max 2140 is calculated using the same values as T.sub.max=M.sub.L+3. The lower right 2170 and upper right 2160 thresholds are calculated using equivalent formulas and the right background area mean value M.sub.R 2165 and its standard deviation . In the example shown in FIG. 21, the values are approximately M.sub.L=110, M.sub.R=90, and =3.3 for both M.sub.L and M.sub.R. The lower and upper band thresholds for all pixels on row 2110 between the two bands 2120 and 2130 are obtained by interpolation. The interpolation between the upper band thresholds 2140 and 2160, and between the lower band thresholds 2150 and 2170, respectively, can be a linear interpolation, as shown in FIG. 21, or another function, such as a quadratic interpolation, if such a function yields a higher resolution binary image.

[0096] 17. The highlighted stippled area 2050 shown in FIG. 20 is the foreground of a first binary image of the grayscale component image shown in FIG. 21, obtained after applying the dynamic thresholding technique described above to all rows of the grayscale component image shown in FIG. 21. The binary image foreground 2050 has a lot of disconnected segments, which sometimes scatter away from the area of the real column. In this case, it is hard to connect those segments to get the column shape. The connections can be made by using a second binary image.

[0097] 18. Create a second binary image by using an isotropic gradient of the same grayscale component image and apply static thresholding based on the gradient background statistic. The isotropic gradient image is assumed to have a zero mean value, so the standard deviation is calculated using only the selected left and right sets of rectangles. A statistical three sigma () rule is used to create the threshold for binarization. The resulting second binary image is shown in FIG. 22 where it is overlaid on top of the isotropic gradient (dark) image.

[0098] 19. Combine the two binary images: the first from the grayscale component image and the second from the isotropic gradient image (above) to produce a composite binary image. The combination enables creating a more complete foreground binary image from disconnected foreground segments. FIG. 23 shows the combined binary image overlaid on the gradient image.

[0099] 20. For further column shape detection and filling of the foreground cavities, two artificial foreground stripes 2410 and 2420 are added, one on the top 2410 and another one on the bottom 2420 (by setting binary pixels to value of 1). FIG. 24 shows the resulting binary overlaid on top of the original grayscale image.

[0100] 21. Morphological operations and a hole filling operation are used to fill foreground holes and smooth rough edges. See Gonzalez, pages 528-536. In FIG. 25, the resulting binary image is shown overlaid on top of the grayscale image.

[0101] 22. The two artificial stripes are removed by setting binary pixels on the top and bottom horizontal edges to background value (0, transparent). Then, a sieving filter is applied for removing small features that account for noise (step 806 in FIG. 8B). The resulting binary image yields a good match between the foreground overlay and the column shape (see FIG. 26), thereby detecting the location of a column of the liquid sample and location of the optical beam path from the composite binary image.

[0102] 23. Connected foreground object(s) are extracted by using a connected component extraction algorithm (see Gonzalez, page 536), thereby evaluating the integrity of the column of the liquid sample from the composite binary image. Normally, just one object matches the normal liquid column. If there are two or more objects, then it is a broken column case (see FIG. 32 for an example). If no objects are found, then it is an empty compartment case. In the latter two cases no further analysis can be done and an error code is reported to notify the operator (steps 807 and 808 in FIG. 8B).

[0103] 24. Calculate the area of the column shape (in pixels) from the detected binary object 2610, as shown in FIG. 26. Additionally, calculate optical path left edge position X.sub.L 2620 and right edge position X.sub.R 2630 by finding the maximum area rectangle that is enclosed into the detected column binary object 2610. The difference X.sub.RX.sub.L is the light path diameter (step 809 in FIG. 8B). The calculated optical beam center X.sub.C 1840, a theoretical maximum optical beam radius R.sub.O (in one embodiment, R.sub.O is equal to about 1/20th of the calculated capital diameter 1110), the column X.sub.L 2620 and X.sub.R 2630 (the enclosed rectangle vertical edge positions), are used for calculating the off-center parameter. The skewedness and off-center conditions of the column of the liquid sample are detected from the composite binary image if the distance min(X.sub.C+R.sub.O, X.sub.R)max(X.sub.CR.sub.O, X.sub.L) is less than 1.5R.sub.O. The theoretical light beam path of radius R.sub.O and centered at X.sub.C should fit into the calculated column light path (rectangle) with the overlap of at least 1.5 light beam radius width.

[0104] Alternative image analysis methods for liquid column shape detection include detecting the object of interest, such as the liquid column, by using edge detection operators for extracting the object contour. These operators are based on computing the difference in pixel intensity between light and dark areas of the grayscale image. There are several basic edge detection (gradient) operators that can be applied: Sobel, Laplacian-of-Gaussian, Roberts, Prewitt, or a composite Canny algorithm. The last consists of several steps, including noise suppression and dynamic thresholding/binarization. See Canny, J., A Computational Approach To Edge Detection, IEEE Trans. Pattern Analysis and Machine Intelligence, Vol. 8(6):679-698, 1986 (hereinafter Canny), the disclosure of which is hereby incorporated by reference in its entirety (however, where anything in the incorporated reference contradicts anything stated in the present application, the present application prevails).

[0105] However, all of these gradient-threshold methods can fall short when dealing with blurry and noisy images, such as an example of a blurry image of a liquid column shown in FIG. 33. The resulting gradient image obtained after applying the Sobel edge detector operator to the blurry image shown in FIG. 33, shown in FIG. 34, reveals the edges (light stripes) of the column but they are washed out on the bottom part. After thresholding/binarization the resulting binary image (stippled overlay) shown in FIG. 35 contains several discontinuous parts with additional spurs. The reason for the failure to detect a contiguous contour (shape) is that the gradient or edge detector operators operate in local image areas (windows), such as 33, 55, etc., pixel wide areas and thus these operators miss high-level structural information about the object contour shape.

[0106] The active contour tracking (also referred to as snakes) method can be used to address the disconnected or noisy (spurious) contour outcome. See Kass, M.; Witkin, A.; Terzopoulos, D. Snakes: Active contour models International Journal of Computer Vision, Vol. 1 (4): 321, 1988, and Chenyang Xu Snakes, Shapes, and Gradient Vector Flow IEEE Transactions on Image Processing, Vol. 7 (3), 1998, the disclosures of which are hereby incorporated by reference in their entirety (however, where anything in the incorporated references contradicts anything stated in the present application, the present application prevails). The active contour tracking method is a combination of the edge detector operators followed by the contour tracking of the result (binarized) that uses the contour curve properties such as its continuity and smoothness. The active contour method is based on the idea of using an energy functional, which represents a weighted combination of internal and external forces that are applied to the contour curve. The internal forces are governed by the physical properties of the contour (elasticity and bending) while the external forces come from the image properties (gradient). The problem is solved by finding an optimum (minimum) of the energy functional. The total energy functional is defined as the following definite integral over the whole range of the parameterized contour v(s)=v(x(s), y(s)), where s belongs to C[0,1]

E*.sub.v=.sub.CE.sub.ac(v(s))ds=.sub.CE.sub.in(v(s))+E.sub.ex(v(s))ds

[0107] where E.sub.in(v(s)) represents the internal energy of the active contour due to elasticity and bending, and E.sub.ex(v(s)) represents the external (image) forces that are applied to the contour. Internal energy is defined as the following two-term sum:

E.sub.in=(|v(s)|.sup.2+|v(s)|.sup.2)/2

[0108] The first-order term, which is controlled by the coefficient, adjusts the elasticity of the active contour. The second-order term, which is controlled by the coefficient, adjusts the stiffness of the active contour. In other words, the first part keeps the active contour short (discourages stretching), while the second part keeps it straight (discourages bending).

[0109] Given a grayscale image L(x,y), which represents a function of intensity in each (x,y)-position of the image, the image (external) force is chosen to lead an active contour toward the object edges and can be represented by two functionals (see Canny):

E.sub.ex.sup.(1)=|L(x,y)|.sup.2

E.sub.ex.sup.(2)=|[G.sub.(x,y)*L(x,y)]|.sup.2

where G.sub.(x, y) is a two-dimensional Gaussian function with standard deviation , is a gradient operator and * denotes the convolution operator. In other words, E.sub.ex.sup.(2) represents a gradient of the smoothed L(x, y) image.

[0110] In the case of the binary image B(x,y) the external forces can be formulated as the following:

E.sub.ex.sup.(1)=.sup.B(x,y)

E.sub.ex.sup.(2)=G.sub.(x,y)*.sup.B(x,y)

where .sup.B(x,y) represents an inverted binary image.

[0111] The parameter controls the smoothness of either the grayscale or the binary imagethe larger the parameter, the blurrier the images and their object edges are. The purpose of the parameter is to extend the search range for the optimization of the active contour.

[0112] The minimum of E*v can be found using the Euler-Lagrange equation:

v(s)v(s)E.sub.ex(v(s))=0

[0113] Let's denote F.sub.in=v(s)v(s) and F.sub.ex=E.sub.ex(v(s)), and then the latter equation can be re-formulated as a force balance equation:

F.sub.in+F.sub.ex=0

[0114] The F.sub.in term represents the internal force that discourages stretching and bending while the external force F.sub.ex pulls the active contour toward the desired image edges. Solving the above equation is accomplished with the gradient descent method, which requires converting the active contour v into a function of time v(s, t). Then the partial derivative of v(s, t) with respect to t can be applied to both sides of the Euler-Lagrange equation. After several iterations when the active contour (snake) has converged to a minimum, its derivative with respect to time becomes zero, and the equation is solved.

v(s,t)/t=v(s,t)v(s,t)E.sub.ex(v(s,t))

[0115] A numerical solution the above equation can be found by discretizing the parameters s and t and solving the equation numerically.

[0116] There is a modification of the active contour model where the curve v is defined implicitly as function of a new parameter r, i.e., s=(r). Then the problem can be reformulated in terms of a geodesic form (GAC), which states that the active contour optimization can be expressed as finding a curve of minimal (geodesic) length under the defined constraints. See Caselles, V., Kimmel, R., and Sapiro, G., Geodesic Active Contours, International Journal of Computer Vision, Vol. 22(1): 61-79, 1997, the disclosure of which is hereby incorporated by reference in its entirety (however, where anything in the incorporated reference contradicts anything stated in the present application, the present application prevails).

[0117] In another aspect, the image analysis techniques described below are concerned with detecting inclusions (in the form of bubbles) and scattered light that can affect the spectral measurements. Two parameters can be measured that are found to be useful in combination: a bubble presence score and scattered light intensity. FIG. 27 shows an example of the image with a significant bubble 2710 near the top of the image.

[0118] An isotropic gradient image of the grayscale component image is used for bubble inclusion detection (see FIG. 28). A special ring detection filter was designed to measure the bubble score. FIGS. 29A and 29B show an example of the ring detection filter with a kernel of 99 pixel size. FIG. 29A shows a 3D contiguous representation of the filter while FIG. 29B shows its digital approximation as a 99 matrix (for display purposes the matrix values are multiplied by 100). The 3D graph shown in FIG. 29A illustrates a raised circular plateau around its perimeter and a negative circular cavity area in the center. The rest of the kernel is set to the intensity value of 0. Its shape correlates with ring-like features on the isotropic gradient image that enables measuring their intensity score, thereby detecting any bubble in the column of the liquid sample.

[0119] A simplified explanation of the ring detection filter is presented in FIGS. 30A-1 through 30B-2, where the ring filter is applied to a single image row (a grey line 3010 near the top of the image). The gradient image with its row profile (the graph shown in FIG. 30A-2) is shown in FIG. 30A-1. The ring filter (top right) applies to all pixels in the row by using the formula shown below the ring filter, where K.sub.j is a kernel coefficient in kernel position j, P.sub.i is a pixel intensity in row position i and F.sub.i is a result value in position i. The process is applied to all image rows. FIG. 30B-1 shows the resulting image with its cross-section row 3020. Note the result peak shown in FIG. 30B-2 that corresponds to the ring-shape feature on the gradient image (marked with an arrow on top).

[0120] The following steps form the ring/bubble presence score calculation algorithm used for detecting any bubble in the column of the liquid sample using both the grayscale component image and the binary image of the light scatter component image:

[0121] 1. apply successive ring detection filtering to a grayscale isotropic gradient image obtained from the grayscale component image (FIG. 28).

[0122] 2. accumulate the result into the cumulative image score.

[0123] 3. use the calculated light path rectangle (2610 obtained in step 24 above, shown in FIG. 26) to limit the ring detection filter application area.

[0124] 4. start from a minimal ring filter size (33) and increase it by 2 (the next is 55) and so on up to a predefined maximum filter diameter (1515, for instance) to cover all possible bubble sizes. While calculating, skip the scores that are below a certain threshold to avoid accumulating values due to noise. The noise threshold is calculated based on statistics of the background rectangular areas for the gradient image (step 810 in FIG. 8B).

[0125] 5. extract the same ROI portion (as that of grayscale component image shown in FIG. 19) from the light scatter component image to estimate the light scatter intensity in the area of the found column shape. FIG. 31 shows a workflow for the light scatter component intensity measurement. Note the noise that is removed and that only bubble scatter is measured.

[0126] 6. as shown in the workflow in FIG. 31, the light scatter component image is binarized to form a binary light scatter component image using a threshold that is calculated using the calculated standard deviation for the selected background rectangular areas (a zero mean is assumed because no blue color is expected to be present in a normal case). Again, a statistical three sigma (a) rule is used to set the binarization threshold. Then a morphological closing (morphological filter) is applied followed by a sieve filter for removing small particles. Then the binary foreground and the calculated light path rectangle are used for limiting the area where the blue component intensity pixels are added up. Then the resulting intensity sum is normalized to make it significant for the value of one and higher (step 811 in FIG. 8B).

[0127] 7. check the bubble presence score and the scattered light intensity score. If both scores are greater that one, then report the defect (steps 812 and 814 in FIG. 8B).

[0128] To distinguish a bubble in the optical beam path from a bubble out of the optical beam path, the calculated optical beam center X.sub.C 1840, a theoretical maximum optical beam radius R.sub.O (in one embodiment, R.sub.O is equal to about 1/20th of the calculated capital diameter 1110), the column X.sub.L 2620 and X.sub.R 2630 (calculated optical path edges) are used. The calculation area for both the bubble detection filter and the scatter light score is limited to the following: left limit is max(X.sub.CR.sub.O, X.sub.L) and the right limit is min(X.sub.C+R.sub.O, X.sub.R). Limiting the calculation area enables performing the bubble and light scatter score measurement only in the part of the image that is known to distort the spectral measurement.

[0129] Images of a liquid column with bubbles, as shown in FIG. 36, suggest that bubble edges (contours) represent circular or elliptical shapes. Note that some of the bubbles can be occluded because the image is a 2D (projection) image of a 3D liquid column with bubbles. For circular object finding, the Hough transform can be used as an alternative method of bubble detection. See Filiberto P., Recognition of Partial Circular Shapes from Segmented Contours, Computer Vision And Image Understanding, Vol. 63(2), 334-342, 1996, the disclosure of which is hereby incorporated by reference in its entirety (however, where anything in the incorporated reference contradicts anything stated in the present application, the present application prevails). Advantages of the transform include its insensitivity to noise and robustness to shape distortion due to occlusion etc.

[0130] Let's consider only circular shapes, although the Hough transform can handle elliptical shapes as well. The Hough transform applies to the binarized (thresholded) image of the gradient image or a binary image of the detected bubble edges. For instance, the Canny operator can be used for edge detection and for thresholding it to a binary image.

[0131] The standard circle equation has the following form:

(xa).sup.2+(yb).sup.2=r.sup.2,

where r is the radius of the circle and (a,b) is a coordinate of the center of the circle.

[0132] The Hough transform applies to the digital form of the circle equation, where all parameters are discrete: x and y are indices of a column and row of a matrix of binary (0 or 1) pixels, parameters a and b are also indices (relative positions) of circle centers and r spans through the possible radii of circles that fit into the image and bound to physical objects of interest (bubbles in this case). The radii usually start from a value greater than one, since an approximation of a radius-one circle on the digital image is too rough (it represents a square). Then, every binary edge (contour) pixel (x.sub.i, y.sub.i) can be transformed into an approximation of a circular cone in the 3D (a, b, r) parameter space. If all contour points lie on a circle, then all its correspondent cones will intersect at a single point (a.sub.i, b.sub.i, r.sub.i) corresponding to the parameters of the circle.

[0133] Since the space is digital, the cones that satisfy the digital form of the circle equation will not intersect at one pixel, but instead represent a small cluster of pixels with a Gaussian-like density distribution, whose center (most dense value) is the resulting (a.sub.i, b.sub.i, r.sub.i) circle triplet. In order to implement the distribution space, an additional voting (integer) value v is needed and the result of the transformation is a 3D matrix of voting values:

V=v(a,b,r),

where a spans though all image columns, b spans through all image rows and r spans through all possible radii of the objects of interest.

[0134] The final and most challenging part of the Hough algorithm is finding the points of local maxima in the resulting matrix V (parametric space). Usually, it requires applying an additional filter for the final matrix V. The resulting points of local maxima can be found by applying a threshold to the filtered matrix V, and they represent all possible circular objects. Because of the voting technique, the algorithm works well even for incomplete or noisy images.

[0135] For the image analysis summary report shown in Table 1, the following parameters with exemplary values are displayed to the operator of the instrument:

TABLE-US-00001 TABLE 1. Image analysis report Off-Center Offset (pixels) 5 Optical Path Diameter (pixels) 208 Light Scatter Intensity Score 21.9 Bubble Presence Score 36.0 Column Feature Area (pixels) 25554 Average ROI Pixel Intensity 148 Column Optical Path Length (pixels) 119

[0136] Off-center Offset: shows the column shape center offset, in pixels (see FIG. 11). Note that pixels are easy to translate to metric units because the physical dimensions of the capital and pedestal are known. The calculated light path rectangle 2610 shown in FIG. 26 and described in step 24 above is used for the off-center offset calculation.

[0137] Optical Path Diameter: the calculated column light/optical path (enclosed cylinder) diameter, in pixels. See step 24 above for details of its calculation.

[0138] Light Scatter Score: measured light scatter normalized intensity, in arbitrary fractional units; a value of 1 and greater usually indicates bubble/inclusions defects. The calculation is shown in step 4 of the ring/bubble presence score calculation algorithm above.

[0139] Bubble Presence Score: in arbitrary units, a value of more than 1 indicates the presence of bubbles. The bubble presence score is used in combination with light scatter score to identify bubbled (defective) columns. For the parameter calculation details, see step 4 of the ring/bubble presence score calculation algorithm above.

[0140] Column Feature Area: measured area of the calculated column shape in pixels.

[0141] Column Optical Path Length: measured height of the calculated light path rectangle in pixels, described in step 24 above.

[0142] Average ROI Pixel Intensity: average image intensity (between 0 and 255), that is useful for detecting underexposed or overexposed grayscale images and adjusting the binarization thresholds.

[0143] The liquid column analysis algorithm produces the following software completion codes:

[0144] Undefined: initial value that means either the analysis was interrupted or failed during initial ROI extraction stage (abnormal condition);

[0145] OK: normal column, expect good spectral reading (step 813 in FIG. 8B);

[0146] Defective column: check the off-center value, the bubble presence score and the light scatter score to identify the reason (step 814 in FIG. 8B);

[0147] Empty Compartment: no liquid column was detected (step 808 in FIG. 8B);

[0148] Broken Column: no connection between the interface surfaces (step 808 in FIG. 8B);

[0149] Column Is Too Short: too short a distance between the interface surface (abnormal condition);

[0150] Poor Background: image background quality is too poor for the analysis (abnormal condition) (step 803 in FIG. 8B).

[0151] While the invention has been described with reference to exemplary embodiments, it will be understood by those skilled in the art that various changes may be made and equivalents may be substituted for elements thereof without departing from the scope of the invention. In addition, many modifications will be appreciated by those skilled in the art to adapt a particular instrument, situation or material to the teachings of the invention without departing from the essential scope thereof. Therefore, it is intended that the invention not be limited to the particular embodiment disclosed as the best mode contemplated for carrying out this invention, but that the invention will include all embodiments falling within the scope of the appended claims.