| Tasks: The purpose of this exercise is to identify a pool of “good” features for the crystallography classification problem. The tasks are 1. Propose at least four features for distinguishing between Class 1 (clear drop) and Class 9 (excellent crystals). 2. Develop an automatic algorithm (executable in Matlab) that extracts these features from a raw image. The algorithm should take care of strong edges and shadows caused by the boundary of the disk. 3. Develop a simple classifier using multivariate Gaussians for the two class-conditional densities, and test the performance of your features by “training” and “testing” on separate halves of the available data. |
| Image Pre-processing | Proposed Features | Decision Classifier | Matlab implementation |
| Finding the dish | Emphasizing the crystals |
One way to detect crystals in an image would be to look for areas that contained high contrast and/or edges. Unfortunately, the dish and the shadows cast by the dish are also sources of high contrast and strong edges. We would like to minimize these sources of contrast and edges before proceeding with our feature extraction. If knowledge of the exact location of the dish in the image was known, then we could simply focus on the interior of the dish when extracting our features. We have developed an algorithm for estimating the location of the dish under the assumption that it is a perfect circle. Our algorithm starts by estimating the edge of the dish at a set of 20 points in the image. This is accomplished by looking at five horizontal and five vertical slices of the image, where each slice results in two estimated points. The slices are located at 6/16, 7/16, 8/16, 9/16, and 10/16 of the length of each axis. The vertical slices contain information pertaining to the top and bottom of the dish, while the horizontal slices contain information for the left and right sides of the dish. As an example of how this information is extracted by our algorithm, consider the third vertical slice, located 8/16 (or 1/2) of the distance across the width (i.e. number of columns) of the image. This slice is a vector with a length equal to the height (i.e. number of rows) of the image. The estimated top (bottom) of the dish in this slice is set equal to the location of the maximum of the absolute value of the finite difference of the first (last) quarter of the slice. This procedure is repeated for the remaining vertical and horizontal slices, which results in 20 estimates of the edge of the dish. Then a circle is fitted, in a least-squares sense, to the points. Since this procedure is looking for the strongest edge within certain regions of the image, it can be thrown off by dark shadows. However, we have found that the least-squares fit tends to average out these sorts of errors.
Below are examples of the resulting dish locations for class 1 (top) and class 9 (bottom) images.
![]() |
![]() |
| back to top | home |
Step 2: Emphasizing the crystals
Once the location of the dish has been estimated, we are free to begin feature extraction, but only using the portion of the image within the dish. But we still must be concerned with shadows which are cast into the dish. These shadows may increase the contrast and/or create an edge in the image which could be picked up by one of our features. This may cause an empty dish to look more like a dish with crystals in it. Also, in the future we may be using more advanced features which attempt to discriminate between great and okay crystals (i.e. 9 vs. 6). Any shadows cast into the dish and over a crystal may mask certain features (e.g. fine 3-D structure) which could be critical to this discrimination. So we would like to further preprocess our images so that only the crystals remain. Since the crystals are mainly defined by edges of a certain size, some sort of algorithm which keeps these edges and suppresses the rest of the image is desired. Our first attempt at such an algorithm consists of a single morphological operation followed by an image subtraction. Specifically, the processed image consists of the difference between the dilated version of the image and the original image. To further enhance the contrast of this processed image, we use the Matlab command 'imadjust'. We have found this process to be very good at suppressing everything in the image except for crystals and the edge of the dish.
Examples of images after pre-processing steps (the first two are the same images shown in step 1):
![]() |
![]() |
![]() |
![]() |
| back to top | home |
| Minimum pixel intensity | Range of pixel intensity | Global sum of edges | Global sum of local variances | Number of objects |
Rationale: When crystals reflect the camera light, they result in extremely bright or dark pixels at the edges, whereas in a clear drop (class 1) the contrast is lower. This suggests that a class-1 image would have lower contrast and so the pixel intensity should fall in the middle of the range. Class-9 images, with sharp edges from the crystals, would have pixel intensity at both ends of the range.
Feature Extraction: This feature is obtained without any pre-processing of the image. First, the histograms of the image intensity are plotted. The histogram bins with less than 5000 pixels (0.35% of total pixels in the image) are then thresholded. Finally the bin with minimum pixel intensity value is selected.


| back to top | home |
Rationale: Again, we use the fact that crystals have sharp edges in the image. Hence, we expect that class-9 images would have a larger range of image intensity than class-1 images.
Feature Extraction: The range of pixel intensity is found by first finding the histogram of the image intensity. Then, the histogram bins with less than 5000 pixels (0.35% of total pixels in the image) are thresholded. The range is the difference between the minimum and maximum pixel intensity (see figure in description of feature 1).

| back to top | home |
Rationale: Class-9 images would have more edges due to the edges of the crystal, while class-1 images would only have edges from the shadows inside the disk. Occasional edges of the disk that does not get removed during the pre-processing step may be present. However, if we sum the total number of pixels detected as edges for each image we should see a larger value for class 9.
Feature Extraction: We use Matlab built-in function ‘edge' to detect edges in the image. The output of the function is a 0-1 logical matrix where 1 is assigned to a pixel belonging to an edge and 0 otherwise. Then we sum all the 1's in the image to get the value that measures how many edges there are in the image.
Example of the detected edge in the image (original image shown in pre-processing step 1):


| back to top | home |
Rationale: The idea is extended from the previous feature that crystals have sharp edges. This is translated into a large value of local variances. Hence, class-9 images which have more edges would have a large local variance. On the other hand, class-1 images would have smooth regions of intensity which would give small local variances.
Feature Extraction: The local variance is calculated for the intensity values at each pixel and its 8 nearest neighbors.
Example of the local variance map in the image (original image shown in pre-processing step 1):


| back to top | home |
Rationale: Besides simple global metrics for contrast, entropy, etc. there are more complicated methods for determining whether or not there are crystals in an image. At a human level, an obvious metric would be the number of crystals in an image. Since we are only using features obtainable by a computer, we will instead attempt to determine how many objects are located within a dish by segmenting the image. The idea is that a dish full of crystals would have lots of segmented objects, while an empty dish, even one with a shadow, would only have a few. A downside of this feature is that a dish with single good crystal in it would probably be classified as empty. More advanced features are needed to screen out these cases.
Feature Extraction: We convert the pre-processed grayscale image to black and white via thresholding. Next, built-in Matlab functions 'bwareaopen' and 'imfill' are used to eliminate clutter and to fill in any small holes. Finally we use 'bwlabel', which segments and labels each distinct object in the image. From this procedure, we obtain the estimated number of objects contained in the dish.
Examples of the segmented and labeled images (same images shown in the pre-processing steps):
![]() |
![]() |
![]() |
![]() |

| back to top | home |
Based upon the performance of our features on the training data (class1: images 0-51, class9: images 0-11), there are no seperations between the two classes in feature 1 (minimum pixel intensity) or feature 2 (range of pixel intensity). Feature 3 (global sum of edges) and feature 4 (global sum of local variances) show some separations, but neither are as good as the separation obtained using the number of objects (feature 5). Therefore, we have chosen to implement a 1-D classifier with Gaussian class-conditional densities having unequal means and variances, equal prior probabilities, and unequal costs (c 91 = 2, c 19 = 5) using only the number of objects as our feature. We have used the training data to estimate the means and variances of class 1 and class 9 for this feature. This gives the decision boundary at approximately 10.8. Specifically, if an image has more than 10 objects it is classified as class 9 (excellent crystals), otherwise it is classified as class 1.

We test our classifier with the remaining images (class1: images 52-104, class9: images 12:24 ). Below is the plot of the test data. Our classifier correctly classifies ALL of the class-9 test data and gives 5 (out of 52) false-alarm errors on the class-1 test data.

| home |
The algorithm for extracting all of the proposed feature is implemented in Matlab. The input of the main function getfeatures is the original image (obtained by calling the function imageconvert). The function pre-process the image as described previously, followed by feature extractions of all the proposed feature. The output of the function is a feature vector of length 5. The m-files are available in the zip file below.
***NOTE*** The m-files provided are written to work with MatLab 7. Error may occur if you run the codes with older versions of MatLab.