From my past blogs, I have discussed morphological operations, blob detection, and connected components. In this blog, we will apply these techniques to extract features in images without the use of Neural Networks.
For our application, the dataset is composed of 27 grayscale images of leaves for five different types of plants. Each image contains a different number of leaves with the same plant type. Our goal is to extract relevant features in classifying the type of plant each leaf belongs to in our image. We will then feed these features to various machine learning models to predict each leaves' plant type and evaluate its accuracy. For simplicity, we will name the five types of plants as plant A, B, C, D, and E. Sample leaves of each plant type are shown in the figure below.
Fig 1. Plant Samples
The general workflow of this application is enumerated below:
Read the image
Binarize the image using a threshold intensity value
Apply morphological operations to clean the image
Detect the connected components in the image
Get the region properties of identified objects
Filter unnecessary objects
Calculate new features per object
Feature selection
Data Preprocessing
Train machine learning model
First, let us import the necessary libraries for our application and define the common function we will use to clean our image.
import os
import re
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from scipy.signal import convolve2d
from skimage.io import imread, imshow
from skimage.color import rgb2gray
from skimage.measure import label, regionprops, regionprops_table
from skimage.transform import rotate
from skimage.morphology import (erosion, dilation, closing, opening,
area_closing, area_opening)
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.svm import LinearSVC
from sklearn.svm import SVC
from sklearn.naive_bayes import GaussianNB
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.ensemble import RandomForestClassifier
from glob import glob
# Constants
plant_prefix = 'Leaves/plant'
sel_h = np.array([[0, 0, 0],
[1, 1, 1],
[0, 0, 0]])
sel_v = np.array([[0, 1, 0],
[0, 1, 0],
[0, 1, 0]])
sel_c = np.array([[0, 1, 0],
[1, 1, 1],
[0, 1, 0]])
blur = (1 / 16.0) * np.array([[1., 2., 1.],
[2., 4., 2.],
[1., 2., 1.]])
def multi_dil(im, num, sel=None):
for i in range(num):
im = dilation(im, sel)
return im
def multi_ero(im, num, sel=None):
for i in range(num):
im = erosion(im, sel)
return im
# Region Properties to extract
properties = ['area', 'centroid', 'convex_area',
'bbox', 'bbox_area', 'eccentricity',
'equivalent_diameter',
'extent', 'filled_area',
'major_axis_length', 'minor_axis_length',
'mean_intensity',
'perimeter', 'orientation', 'solidity']Now that we have imported, defined, and assigned the necessary libraries, functions, and variables, we can now segment each leaf in the dataset. We binarize the images using a pre-defined threshold to remove unnecessary objects in our image. However, some leaves are close to each other and appear as one leaf. Morphological operations can be used to solve this issue and separate conjoined leaves. After isolating each leaves with each other, we use the label and regionprops_table function from the skimage library. The regionprops_table function extracts each connected components' properties in our image, such as area, centroid, mean intensity, etc. However, not all the objects detected are leaves. We used the convex_area of each object to filter out objects that did not meet specific standards to be classified as a leaf. We can also use the orientation property to derive new features not available in the regionprops_table function like the circularity and length of each leaf to increase the accuracy of the machine learning models. The code for segmenting and extracting features for plant a is shown in the figure below.
plant_a_files = glob(f'{plant_prefix}A*.jpg')
clean_plant_a = [] # List of two tuples (img, regionprops)
fig = plt.figure(figsize=(20, 20))
file_count = len(plant_a_files)
thres = 0.4
no_col = 6
no_row = int(np.ceil(file_count * 2 / no_col))
gs = gridspec.GridSpec(no_row, no_col)
for i, file in enumerate(plant_a_files):
img = imread(file)
fn = file.split('\\')[-1].split('.')[0]
ax0 = fig.add_subplot(gs[i * 2])
ax1 = fig.add_subplot(gs[i * 2 + 1])
ax0.axis('off')
ax1.axis('off')
# Display Threshold Image
ax0.imshow(img)
ax0.set_title(fn)
# Threshold image
img = rgb2gray(img)
img_bw = img < thres
# Morph Image
img_morph = area_closing(area_opening(img_bw, 200), 200)
img_morph = multi_ero(multi_ero(img_morph, 22, sel_v), 2, sel_h)
img_morph = multi_dil(img_morph, 20, sel_v)
# Get region properties of image
img_label = label(img_morph)
df_regions = pd.DataFrame(regionprops_table(img_label, img,
properties=properties))
# Filter regions using area
area_thres = df_regions.convex_area.mean()/3
df_regions = df_regions.loc[df_regions.convex_area > area_thres]
mask_equal_height = ((df_regions['bbox-2'] - df_regions['bbox-0'])
!= img_label.shape[0])
mask_equal_width = ((df_regions['bbox-3'] - df_regions['bbox-1'])
!= img_label.shape[1])
df_regions = df_regions.loc[mask_equal_height & mask_equal_width]
# Compute for Derive features
y0, x0 = df_regions['centroid-0'], df_regions['centroid-1']
orientation = df_regions.orientation
x1 = (x0 + np.cos(df_regions.orientation) * 0.5
* df_regions.minor_axis_length)
y1 = (y0 - np.sin(df_regions.orientation) * 0.5
* df_regions.minor_axis_length)
x2 = (x0 - np.sin(df_regions.orientation) * 0.5
* df_regions.major_axis_length)
y2 = (y0 - np.cos(df_regions.orientation) * 0.5
* df_regions.major_axis_length)
df_regions['major_length'] = np.sqrt((x2 - x0)**2 + (y2 - y0)**2)
df_regions['minor_length'] = np.sqrt((x1 - x0)**2 + (y1 - y0)**2)
df_regions['circularity'] = (4 * np.pi * df_regions.filled_area
/ (df_regions.perimeter ** 2))
# Display segmented image
ax1.imshow(img_label)
ax1.set_title(f'{fn} segmented: {df_regions.shape[0]}')
df_regions['target']='plant_a'
clean_plant_a.append(df_regions)
Fig 3. Plant A segmentation
In segmenting other plant types, we can reuse the steps for segmenting plant A. Now that we have extracted the features from different plant types, feature selection is performed to efficiently train the machine learning model and avoid storing unnecessary data. In building our dataset, we will use the following features:
area - Area of the detected object
convex_area - Area of the detected convex region
bbox_area - Number of bounding box pixels.
eccentricity - Ratio of the focal distance over the major axis length.
equivalent_diameter - Diameter of a circle with the same area as the region.
extent - Ratio of pixels in the region to pixels in the total bounding box
filled_area - Number of pixels of the region will all the holes filled in.
perimeter - Length of the path that surrounds the detected object.
mean_intensity - Average intensity in the detected region.
solidity - Ratio of pixels in the region to pixels of the convex hull.
major_length - Major length of the object.
minor_length - Minor length of the object
circularity - Degree of roundness of the detected object.
To know more about the region properties, check this link.
df_plant = pd.concat([*clean_plant_a, *clean_plant_b, *clean_plant_c,
*clean_plant_d, *clean_plant_e])
req_feat = ['area', 'convex_area', 'bbox_area', 'eccentricity',
'equivalent_diameter', 'extent', 'filled_area',
'perimeter', 'mean_intensity',
'solidity', 'major_length', 'minor_length',
'circularity', 'target']
df_plant = df_plant[req_feat]
df_plant
Fig 4. Extracted Features
Since the features we have extracted have different scales, the dataset is normalized using the Min-Max normalization method. To validate the model's robustness and generalizability, we have also split the dataset into two, training and test dataset.
X = df_plant.drop('target', axis=1)
y = df_plant.target
# Normalize the data
scaler = MinMaxScaler()
X_scaled = scaler.fit_transform(X)
# Split the data
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y,
test_size=0.25,
random_state=0)Now that we have extracted and preprocessed our data, we are ready to train our models. Different machine models are used to evaluate which model performed best given the features we have extracted in the dataset. The top-predictors for each plant type are also monitored to determine which features are significant in classifying the plant type.
#Training
startTime = time.time()
knn = KNeighborsClassifier()
knn.fit(X_train, y_train)
elapsed_time_knn = time.time() - startTime
startTime = time.time()
LR2 = LogisticRegression(penalty='l2', max_iter=10000, n_jobs=-1)
LR2.fit(X_train, y_train)
elapsed_time_LR2 = time.time() - startTime
startTime = time.time()
LR1 = LogisticRegression(penalty='l1', max_iter=10000, solver='liblinear')
LR1.fit(X_train, y_train)
elapsed_time_LR1 = time.time() - startTime
startTime = time.time()
LSVC2 = LinearSVC(penalty='l1', dual=False, max_iter=10000)
LSVC2.fit(X_train, y_train)
elapsed_time_LSVC2 = time.time() - startTime
startTime = time.time()
LSVC1 = LinearSVC(penalty='l2', max_iter=10000)
LSVC1.fit(X_train, y_train)
elapsed_time_LSVC1 = time.time() - startTime
startTime = time.time()
NSVC = SVC(kernel='rbf')
NSVC.fit(X_train, y_train)
elapsed_time_NSVC = time.time() - startTime
startTime = time.time()
Bayes_Gaussian = GaussianNB()
Bayes_Gaussian.fit(X_train, y_train)
elapsed_time_Bayes_Gaussian = time.time() - startTime
startTime = time.time()
DT = DecisionTreeClassifier(max_depth=5)
DT.fit(X_train, y_train)
elapsed_time_DT = time.time() - startTime
startTime = time.time()
RF = RandomForestClassifier(max_depth=5, n_estimators=100)
RF.fit(X_train, y_train)
elapsed_time_RF = time.time() - startTime
startTime = time.time()
GB = GradientBoostingClassifier(max_depth=5, n_estimators=100)
GB.fit(X_train, y_train)
elapsed_time_GB = time.time() - startTime
#Prediction
startTime = time.time()
y_pred_knn = knn.predict(X_test)
test_time_knn = time.time() - startTime
startTime = time.time()
y_pred_LR2 = LR2.predict(X_test)
test_time_LR2 = time.time() - startTime
startTime = time.time()
y_pred_LR1 = LR1.predict(X_test)
test_time_LR1 = time.time() - startTime
startTime = time.time()
y_pred_LSVC2 = LSVC2.predict(X_test)
test_time_LSVC2 = time.time() - startTime
startTime = time.time()
y_pred_LSVC1 = LSVC1.predict(X_test)
test_time_LSVC1 = time.time() - startTime
startTime = time.time()
y_pred_NSVC = NSVC.predict(X_test)
test_time_NSVC = time.time() - startTime
startTime = time.time()
y_pred_Bayes_Gaussian= Bayes_Gaussian.predict(X_test)
test_time_Bayes_Gaussian = time.time() - startTime
startTime = time.time()
y_pred_DT= DT.predict(X_test)
test_time_DT = time.time() - startTime
startTime = time.time()
y_pred_RF= RF.predict(X_test)
test_time_RF = time.time() - startTime
startTime = time.time()
y_pred_GB= GB.predict(X_test)
test_time_GB = time.time() - startTime
cols = ['Model','Train Accuracy',
'Test Accuracy', 'Training Time', 'Testing Time',
"Plant A Top Predictor",
"Plant B Top Predictor",
"Plant C Top Predictor",
"Plant D Top Predictor",
"Plant E Top Predictor",
"Top Predictor Importance"]
df_result = pd.DataFrame(columns=cols)
df_result.loc[0] = ['kNN', knn.score(X_train, y_train),
knn.score(X_test, y_test),
elapsed_time_knn,
test_time_knn,
"", "", "", "", "", ""]
df_result.loc[1] = ['Logistic Regression (L2)', LR2.score(X_train, y_train),
LR2.score(X_test, y_test),
elapsed_time_LR2,
test_time_LR2,
X.columns[np.argmax(LR2.coef_[0])],
X.columns[np.argmax(LR2.coef_[1])],
X.columns[np.argmax(LR2.coef_[2])],
X.columns[np.argmax(LR2.coef_[3])],
X.columns[np.argmax(LR2.coef_[4])],
np.max(LR2.coef_)/np.sum(LR2.coef_)]
df_result.loc[2] = ['Logistic Regression (L1)', LR1.score(X_train, y_train),
LR1.score(X_test, y_test),
elapsed_time_LR1,
test_time_LR1,
X.columns[np.argmax(LR1.coef_[0])],
X.columns[np.argmax(LR1.coef_[1])],
X.columns[np.argmax(LR1.coef_[2])],
X.columns[np.argmax(LR1.coef_[3])],
X.columns[np.argmax(LR1.coef_[4])],
np.max(LR1.coef_)/np.sum(LR1.coef_)]
df_result.loc[3] = ['Linear SVC (L2)', LSVC2.score(X_train, y_train),
LSVC2.score(X_test, y_test),
elapsed_time_LSVC2,
test_time_LSVC2,
X.columns[np.argmax(LSVC2.coef_[0])],
X.columns[np.argmax(LSVC2.coef_[1])],
X.columns[np.argmax(LSVC2.coef_[2])],
X.columns[np.argmax(LSVC2.coef_[3])],
X.columns[np.argmax(LSVC2.coef_[4])],
np.max(LSVC2.coef_)/np.sum(LSVC2.coef_)]
df_result.loc[4] = ['Linear SVC (L1)', LSVC1.score(X_train, y_train),
LSVC1.score(X_test, y_test),
elapsed_time_LSVC1,
test_time_LSVC1,
X.columns[np.argmax(LSVC1.coef_[0])],
X.columns[np.argmax(LSVC1.coef_[1])],
X.columns[np.argmax(LSVC1.coef_[2])],
X.columns[np.argmax(LSVC1.coef_[3])],
X.columns[np.argmax(LSVC1.coef_[4])],
np.max(LSVC1.coef_)/np.sum(LSVC1.coef_)]
df_result.loc[5] = ['Non-linear SVC', NSVC.score(X_train, y_train),
NSVC.score(X_test, y_test),
elapsed_time_NSVC,
test_time_NSVC,
"", "", "", "", "", ""]
df_result.loc[6] = ['Naive Bayes Gaussian', Bayes_Gaussian.score(X_train, y_train),
Bayes_Gaussian.score(X_test, y_test),
elapsed_time_Bayes_Gaussian,
test_time_Bayes_Gaussian,
"", "", "", "", "", ""]
df_result.loc[7] = ['Decision Tree', DT.score(X_train, y_train),
DT.score(X_test, y_test),
elapsed_time_DT,
test_time_DT,
X.columns[np.argmax(DT.feature_importances_)],
X.columns[np.argmax(DT.feature_importances_)],
X.columns[np.argmax(DT.feature_importances_)],
X.columns[np.argmax(DT.feature_importances_)],
X.columns[np.argmax(DT.feature_importances_)],
np.max(DT.feature_importances_)]
df_result.loc[8] = ['Random Forest Classifier', RF.score(X_train, y_train),
RF.score(X_test, y_test),
elapsed_time_RF,
test_time_RF,
X.columns[np.argmax(RF.feature_importances_)],
X.columns[np.argmax(RF.feature_importances_)],
X.columns[np.argmax(RF.feature_importances_)],
X.columns[np.argmax(RF.feature_importances_)],
X.columns[np.argmax(RF.feature_importances_)],
np.max(RF.feature_importances_)]
df_result.loc[9] = ['Gradient Boosting Classifier', GB.score(X_train, y_train),
GB.score(X_test, y_test),
elapsed_time_GB,
test_time_GB,
X.columns[np.argmax(GB.feature_importances_)],
X.columns[np.argmax(GB.feature_importances_)],
X.columns[np.argmax(GB.feature_importances_)],
X.columns[np.argmax(GB.feature_importances_)],
X.columns[np.argmax(GB.feature_importances_)],
np.max(GB.feature_importances_)]
df_result
Fig 5. Accuracy and Top predictors of each model.
In summary, training the plant classifier using the extracted features using traditional image processing techniques in the image, we achieved a validation accuracy of 96.9% using the Gradient Boosting classifier. Using each segmented leaf's regional properties, we can extract the rich features like area, convex_area, bbox_area, eccentricity, equivalent_diameter, extent, filled_area, mean_intensity, perimeter, and solidity.
Comments