On the classification of simple and complex biological images using Krawtchouk moments and Generalized pseudo-Zernike moments: a case study with fly wing images and breast cancer mammograms

Problem 1: fly wing venation patterns for species identification

Generally, identifying a biological specimen down to the species level with certainty requires a certain level of taxonomic expertise. The taxonomist examines morphological characteristics of the specimen physically, and applies expert judgement to classify the specimen. This process is often slow and expensive. Additionally, taxonomists may also be increasingly hard to find in the future, as the number of permanent positions stagnate or shrink as a consequence of lack of funding and training at the tertiary level (Britz, Hundsdörfer & Fritz, 2020).

Traditional morphometric analysis (Marcus, 1990), which mainly captures size variation, or geometric morphometric analysis (Bookstein, 1991; Adams, Rohlf & Slice, 2013), which captures shape variation, are possible quantitative methods that potentially allow a data-driven approach to species identification. Landmark-based geometric morphometrics relies on using homologous landmarks, which can be unambiguously identified on an image. However, depending on the organism of interest, it is possible that few or no homologous landmarks are available, despite the fact that biological shape variation is apparent (e.g. cellular shape, claw shape) to the human observer.

Species identification by analysis of wing venation patterns often leads to correct identification at the species or even subpopulation level because the main source of variation in wing venation patterns is evolutionary divergence between taxa, with only rare and incomplete secondary convergence (Perrard et al., 2014). Currently, there is persistent interest in applying geometric morphometric analysis of wing venation patterns as a basis for identifying forensically important flies (e.g. Sontigun et al. (2017), Sontigun et al. (2019)). Recently, Khang et al. (2021) provided proof-of-concept that the species identities of forensically important flies predicted using wing venation geometric morphometric data with random forests are highly concordant with those inferred from DNA sequence data. Since analysis of whole wing image is likely to yield higher resolution data, we hypothesize that this will yield improved species prediction performance compared to using geometric morphometric landmark data which is relatively low resolution. Indeed, Macleod, Hall & Wardhana (2018) reported encouraging results from an image analysis of fly wing venation patterns using pixel brightness as features. However, their method requires the use of undamaged wings and a standardized protocol to minimize imaging artefacts arising from slide preparation (e.g. bubbles, lighting variation). Therefore, capturing information in image pixel data as translation, rotation and scale-invariant features may improve usability of images without the need to apply a rigid imaging protocol.

Here, we do not consider comparison against Elliptic Fourier Analysis (Kuhl & Giardina, 1982), another shape analysis method, since it is used for shapes that are closed contours, which wing venation patterns are not.

Problem 2: breast mass classification

The classification of breast masses based solely on digital mammograms is a challenging problem, owing to the heterogeneous morphology of breast masses (Aleskandarany et al., 2018). Several researchers who used orthogonal moments to construct image features for the classification of benign and malignant masses reported encouraging findings (Tahmasbi, Saki & Shokouhi, 2011; Narváez & Romero, 2012). Current state-of-the-art deep learning approach to image analysis of breast cancer mammograms gave highly optimistic results. Shen et al. (2019) reported sensitivity of 86.7%, and specificity of 96.1% in the classification of benign and malignant breast masses, using 2478 images in the CBIS-DDSM database (training set size = 1903; validation set size = 199; test set size =376; (Clark et al., 2013; Lee et al., 2017)). Nevertheless, the opacity and plasticity of powerful black-box methods such as deep learning pose challenges to their formal adoption in medical practice (Nicholson Price, 2018). In the end, combining the complementary strengths of features derived from human expert judgement and those from statistical learning models seems to be the most convincing approach (Gennatas et al., 2020).

Here, we hypothesize that integrating the result of statistical learning outcome from image analysis using orthogonal moments with relevant expert features may ameliorate performance deficiencies based solely on image analysis.

Fly wing images

For the problem of fly species identification using images of wing venation patterns, we used species from two forensically important fly families: Sarcophagidae and Calliphoridae (Amendt et al., 2011). Images of wings of male specimens that are of sufficiently good quality for image analysis, and their associated geometric morphometric data from 19 landmarks (Fig. 1) were taken from Khang et al. (2020).

The samples were taken from male flies from the Calliphoridae family (seven species), namely Chrysomya megacephala (Fabricius, 1794) (n = 5), Chrysomya nigripes Aubertin, 1932 (n = 5), Chrysomya pinguis (Walker, 1858) (n = 5), Chrysomya rufifacies (Macquart, 1842) (n = 5), Chrysomya villeneuvi Patton, 1922 (n = 5), Lucilia cuprina (Wiedemann, 1830) (n = 4) and Lucilia porphyrina (Walker, 1856) (n = 5). From the Sarcophagidae family (eight species), we had Boettcherisca javanica Lopes, 1961 (n = 5), Boettcherisca karnyi (Hardy, 1927) (n = 5), Boettcherisca peregrina (Robineau-Desvoidy, 1830) (n = 5), Sarcophaga ruficornis (Fabricius, 1794) (n = 5), Sarcophaga dux Thompson, 1869 (n = 5), Parasarcophaga albiceps (Meigen, 1826) (n = 5), Parasarcophaga misera (Walker, 1849) (n = 5) and Sarcophaga princeps Wiedemann, 1830 (n=5). In total, 74 specimens from 15 different species were used.

Breast cancer images and associated expert features

The DDSM (Digital Database for Screening Mammography) database (Heath et al., 1998, 2000) is a public repository of breast cancer mammograms. Although this database contains a large collection of 2,620 scanned film mammograms, the quality of annotations in the images varies. Examples of errors include wrongly annotated images and lesion outlines that do not form precise mass boundary (Lee et al., 2017; Song et al., 2009). Inclusion of such poor quality images in the training phase of a statistical learning model can weaken model generalizability.

To overcome this problem, a curated subset of images in the DDSM database, known as the CBIS-DDSM (Curated Breast Imaging subset of DDSM) collection (Clark et al., 2013; Lee et al., 2017) was created. Images in this database consist of selected mammograms that have been segmented using an automated segmentation algorithm. The segmented images were evaluated by comparing outlines of mass lesion images with hand-drawn outlines made by a trained radiologist. The CBIS-DDSM collection comprises scanned filmed mammography from 1,566 participants. A patient could have more than one type of lesions (e.g. mass, calcification) in a single mammogram. The images were decompressed and converted to the DICOM format containing updated region of interest (ROI) segmentation, bounding boxes, and pathologic diagnosis. Moreover, a total of 3,568 focused images are also available in this database to cater to studies that do not require the use of full mammogram images but only a focused region of abnormalities. We used 1,151 images from the CBIS-DDSM collection, of which 48% (556/1,151) are of benign class, and 52% (595/1,151) are of malignant class. The images were downloaded from the CBIS-DDSM database in the PNG file format.

Each mammogram image is further associated with five expert features: (i) BI-RADS assessment; (ii) mass shape; (iii) mass margin; (iv) breast density; (v) subtlety rating. The Breast Imaging-Reporting and Data System (BI-RADS) provides a standard for reporting breast examination results based on mammography, ultrasound, and magnetic resonance imaging data. First published in 1992 (American College of Radiology, 1992), BI-RADS has become a standard communication tool for mammography reports globally (Balleyguier et al., 2007), and is now in its fifth edition (Sickles et al., 2013). By standardizing the reporting of mammography results, BI-RADS facilitates communication among radiologists and clinicians, and aids the training and education of junior radiologists in developing countries (Lehman et al., 2001).

There are seven categories in the BI-RADS assessment that can be assigned based on evaluation of the lexicon descriptors or the biopsy findings of a lesion (Sickles et al., 2013). Category 0 indicates that materials are insufficient for evaluation and additional imaging evaluation or prior mammograms for comparison are required. Category 1 is given when no anomalies are found. Category 2 is given when there is evidence of benign tumors such as skin calcifications, metallic foreign bodies, fat-containing lesions and involuting calcified fibroadenomas. If radiologists are unsure of the lesion categorization, a BI-RADS category of 3 is given, and a follow-up over a certain interval of time is done to determine stability of the lesion. The risk of malignancy in this category is considered to be at most 2%. Category 4 is assigned when malignant tumors are suspected. Three subcategories are possible: a, b, and c. The subcategory (a) reflects a subjective probability of malignancy between 2% to 10%. Subcategory (b) reflects a subjective probability of 10% to 50% for malignancy, while subcategory (c) reflects a subjective probability of malignancy ranging from 50% to 95%. Category 5 reflects a strong belief (probability of 95% or more) that a lesion is malignant. When a lesion placed in this category is contradicted by a benign biopsy report, a surgical consultation may still be advised. Finally, a BI-RADS category of 6 is given when a lesion receives confirmation of malignancy from the biopsy result.

A breast mass has two important aspects: shape and margin. BI-RADS lexical descriptors of mass shape include oval, irregular, lobulated, etc. The mass margin describes the shape of the edges of a mass, such as being circumscribed, ill-defined, or spiculated. Breast density (Sickles et al., 2013) is a categorical variable with four categories. Category 1 describes breast composition that is almost entirely fatty. Category 2 indicates the presence of scattered fibroglandular densities. Category 3 indicates breast that is heterogeneously dense. Category 4 indicates breast that is extremely dense, which lowers sensitivity of mammography. Finally, the subtlety rating, which is not part of the BI-RADS standards, is an ordinal variable on a scale of 1 to 5 representing the difficulty in viewing the abnormality in a mammogram (Lee et al., 2017). The scale ranges from 1 for “subtle” to 5 for “obvious”.

Data processing

Feature extraction

Fly wing images

To use geometric morphometric data from the wing images, raw coordinate data from the 19 landmarks were processed using Generalized Procrustes Analysis in the geomorph R package (Version 3.3.2; Adams & Otarola-Castillo, 2013) to produce the translation, rotation, and scale-invariant Procrustes coordinates. To remove the effect of allometry, we used the residuals produced from linear regression of the Procrustes coordinates against the logarithm (base 10) of centroid size (Sidlauskas, Mol & Vari, 2011; Klingenberg, 2016). To remove correlation between the Procrustes coordinates, we applied R-mode principal component analysis (PCA), and kept the first 15 principal components that cumulatively explain 98.7% of the total variation in the data.

For image analysis, since the fly wing images are rectangular, we used higher order moments to ensure that the reconstruction captured image details distal from the image centroid (Fig. 2). We found that Krawtchouk moment invariants (see Appendix) of order 200 was appropriate for extracting image features from the binary images of the fly wings. Thus, 40,000 moment invariant features were obtained. Since the range of values for these features was generally large, we scaled these features using the Z-score, and then applied Q-mode PCA to reduce the dimension of the feature space. The first 60 principal components accounting for 92.6% of total variance were used as features for downstream statistical learning work.

Study design and analysis

Performance evaluation

We used standard metrics for the evaluation of classifier performance. Accuracy is defined as the probability that the predicted class is the same as the true class. For multi-class prediction in the fly images, we did not use sensitivity or specificity, as contextually no particular species is of any special interest. Hence, sensitivity and specificity were considered only in the breast cancer classification. There, sensitivity is defined as the probability of predicting the malignant class, given that a sample is malignant. Specificity is defined as the probability of predicting the benign class, given that a sample is benign.

In the case of the fly images, we used the Bayesian posterior mean of accuracy with uniform prior (see Appendix), and reported the 95% Bayesian credible interval (Brown, Cai & DasGupta, 2001). For the breast cancer images, we reported the mean of accuracy, sensitivity, and specificity from the 10 random instances, along with their associated standard error estimate (sample standard deviation /√10).

Software and computation

For image analysis, we used R version 3.6.1 (R Core Team, 2018) to perform the computations and run the IM R package (Rajwa et al., 2013). Data processing and analyses of the CBIS-DDSM samples and fly wing images were done using a 22 CPU core, 23 GB RAM server running on Ubuntu 16.04.4 LTS at the Data Intensive Computing Centre, Universiti Malaya, Malaysia. For classification using decision trees with kernel model and random forests, we used the GUIDE program (Version 35.2; Loh (2009, 2014)).

Quality of moment feature representation of images

For the fly wings, images reconstructed from KM of order 200 approximated the binary images very well (Fig. 2C). Similarly, images reconstructed from KM (order 163) and GPZM (α = 0) for the breast images also approximated the ROI well. An example showing a malignant mass is given in Fig. 3. For additional examples, see Supplemental File 2.

Classification of fly species

Krawtchouk moment invariants of order 200 captured important patterns of variation in wing venation that led to 94.7% correct classification of 74 training samples (95% Bayesian credible interval for accuracy = [88.8–98.5%]). Indeed, within species variation was substantially smaller for 12 species, compared to between species variation (Fig. 4). In comparison, the use of geometric morphometric data produced classification accuracy of 64.5% (95% Bayesian credible interval for accuracy = [53.5–74.8%]). Baseline prediction accuracy using the majority class was 7.9%. The classification matrices for both cases are given in Supplemental File 3.

Under the random forests model, predicting species using Krawtchouk moment invariants produced 100% accurate classification result, whereas two errors were made when geometric morphometric data was used. Out-of-bag error using Krawtchouk moment invariants was 0% (0/74), compared to 39.2% using geometric morphometric data (29/74).

Classification of benign and malignant breast masses

The mean classification accuracy based solely on image data was about 57% ± 1%, with mean sensitivity about 70% ± 1%, and mean specificity of 43% ± 2%. Baseline prediction accuracy using majority class was 52%. Figure 5 shows, for a particular training-testing instance, the estimated bivariate Gaussian kernel densities in the space of the first linear discriminants derived from KM and GPZM, with test samples superimposed on the plot.

When expert features were used together with P(mal) as a feature, we observed substantial increase of mean accuracy (from 57% to between 68% and 75%; Table 1), mean sensitivity (from 70% to between 80% and 97%), and mean specificity (from 43% to between 51% and 61%, excluding Model IV). Using the set of expert features without P(mal) (Model II) produced the best mean accuracy (75% ± 1%) and best mean sensitivity (97% ± 0%).

Removal of the BI-RADS assessment (Model IV) affected mean specificity substantially, causing it to drop from 51% (Model II) to 40%. If BI-RADS assessment was replaced by P(mal), mean specificity improved slightly, at the expense of mean sensitivity, which dropped from 97% (Model II) to 80%, thus causing a concurrent drop in mean accuracy from 75% to 68%. This suggests that the information in BI-RADS assessment and P(mal) are non-redundant and may complement each other. Finally, these two features interact to produce a model (Model I) where the difference in sensitivity and specificity is minimized (24%) compared to other models (Model II: 46%; Model III: 34%; Model IV: 55%). Thus, the inclusion of P(mal) as a feature seemed to be important for reducing the discrepancy between model sensitivity and model specificity by decreasing the former but increasing the latter.

Table 2 shows summary statistics of feature importance (mean ± standard deviation) for each of the four models. Across all four models, BI-RADS assessment (where used) and mass margin were consistently ranked as the two most important features, and the order of importance (in decreasing importance) for mass shape, subtlety, and breast density was also consistent. Where P(mal) was used in Model I, on average it ranked third in order of importance; in Model III, on average it ranked second, after mass margin. This suggests that P(mal), which summarizes abstract information from image data, provides useful information to the statistical learning model when used together with the expert features.

Orthogonal moments

In this section, we provide sufficient mathematical background for the appreciation of the use of orthogonal moments as feature extractors of images. For further advanced details, we refer readers to Szegö (1975) for the theory of orthgonal polynomials, Yap, Paramesran & Ong (2003) for Krawtchouk moments, Xia et al. (2007) for generalized pseudo-Zernike moments, and Shu, Luo & Coatrieux (2007) for a general introduction to using orthogonal moments for image analysis.

Mathematical preliminaries

Definition 1. (Szegö, 1975) Let p_n(x) be a polynomial in x of order n. For the interval a ≤ x ≤ b, if w(x) is a weight function in x, and δ_nm is the Kronecker delta which is equal to 1 when n = m, and 0 when n ≠ m, then p_n(x) is said to be an orthogonal polynomial associated with the weight function w(x) if it satisfies the condition

$${\int }_a^b{p}_n(x)p_m(x)w(x)dx={\delta }_{nm}.$$

Definition 2. (Oberhettinger, 1964) The hypergeometric function, ₂F₁(a,b;c;z) is a special function defined as the power series

$${}_2{F}_1(a,b;c;z)=\sum _{k=0}^{\mathrm{\infty }}\frac{(a{)}_k(b{)}_k}{(c{)}_k}\frac{z^k}{k!},$$where a,b,c are real numbers, with |z| < 1. The notation (a)_k denotes the Pochhammer symbol for the rising factorial (a)_k = a(a + 1)(a + 2)…(a + k − 1), with (a)₀ = 1.

Krawtchouk polynomials and moments

The Krawtchouk polynomials (Krawtchouk, 1929a, 1929b) are discrete orthogonal polynomials associated with a binomial probability weight function. The Krawtchouk polynomial of order n is denoted by $k_n(x;p,N-1)$, and can be conveniently representing as a hypergeometric function

$$k_n(x;p,N-1){=}_2{F}_1(-n,-x;-N+1;\frac{1}{p})$$where n, x = 0, 1, 2, … , N − 1, N >1, 0 < p < 1. The weighted Krawtchouk polynomial (Yap, Paramesran & Ong, 2003) ${\overline{k}}_n(x;p,N-1)$ is given by

$${\overline{k}}_n(x;p,N-1)=k_n(x;p,N-1)\sqrt{\frac{\omega (x;p,N-1)}{\rho (x;p,N-1)}},$$where $\omega (x;p,N-1)$ is the binomial probability mass function

$$\omega (x;p,N-1)=(\genfrac{}{}{0ex}{}{N-1}{x})p^x(1-p{)}^{N-1-x},$$and $(1-p{)}^{N-1}/\rho (x;p,N-1)=\omega (n;p,N-1)$. For brevity, we will write ${\overline{k}}_n(x;p,N-1)$ as ${\overline{k}}_n(x)$.

Let the Krawtchouk moments matrix of order l be an l × l square matrix, Q. The (n,m) element of Q, denoted as Q_nm, is related to the image intensity function f(x,y) on the two-dimensional discrete domain through

$$Q_{nm}=\sum _{x=0}^{N-1}\sum _{y=0}^{M-1}{\overline{k}}_n\left(x\right){\overline{k}}_m\left(y\right)f(x,y).$$

Consider an image of dimension N × M. If we denote K₁ and K₂ as the l × N and l × M matrix of weighted Krawtchouk polynomials, respectively, with A as the N × M matrix of image intensity functions, then

$$\mathit{Q}={\mathit{K}}_1\mathit{A}{\mathit{K}}_2^T.$$

The orthogonality property of the weighted Krawtchouk polynomials

$${\displaystyle \sum _{x=0}^{N-1}{\overline{k}}_n}\left(x\right){\overline{k}}_m\left(x\right)={\delta }_{nm},$$implies that the product of an l × N matrix of weighted Krawtchouk polynomials with its transpose is the l × l identity matrix. Therefore, the image intensity functions can be reconstructed from $\mathit{Q}$ as

$$\mathit{A}={\mathit{K}}_{\mathbf{1}}^T\mathit{Q}{\mathit{K}}_{\mathbf{2}},$$that is,

$$f(x,y)=\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}{\overline{k}}_n\left(x\right){\overline{k}}_m\left(y\right)Q_{nm}.$$

Generalized pseudo-Zernike polynomials and moments

The pseudo-Zernike polynomials (Bhatia & Wolf, 1954) are polynomials in two variables that form a complete orthogonal set for the interior of the unit circle. The pseudo-Zernike polynomials of order n and repetition m is denoted by V_nm(r,θ), and expressed in polar coordinate form as

$$V_{nm}(r,\theta )=R_{nm}(r)e^{im\theta },$$where i is the complex number $\sqrt{-1}$, and R_nm(r) is the radial polynomial defined as

$$R_{nm}(r)=\sum _{j=0}^{n-|m|}\frac{(-1{)}^j(2n+1-j)!}{j!(n-|m|-j)!(n+|m|+1-j)!}{r}^{n-j},$$where $n=0,1,2,\dots ,\mathrm{\infty },$ and |m| ≤ n. Xia et al. (2007) proposed a generalization of R_nm(r)

$$R_{nm}^{\alpha }(r)=\frac{(n+|m|+1)!}{(\alpha +1{)}_{n+|m|+1}}\sum _{j=0}^{n-|m|}(-1{)}^j\frac{(\alpha +1{)}_{2n+1-j}}{j!(n-|m|-j)!(n+|m|+1-j)!}{r}^{n-j},$$where α > − 1, with R⁰_nm(r) = R_nm(r). The weighted generalized radial polynomial is given by

$${\overline{R}}_{nm}^{\alpha }(r)=R_{nm}^{\alpha }(r)\sqrt{\frac{(2n+\alpha +2)(\alpha +1+n-|m|{)}_{2|m|+1})}{2\pi (n-|m|+1{)}_{2|m|+1}}}(1-r{)}^{\alpha /2},$$leading to the generalized pseudo-Zernike polynomials

$${\overline{V}}_{nm}^{\alpha }(r,\theta )={\overline{R}}_{nm}^{\alpha }(r)e^{im\theta }.$$

The generalized pseudo-Zernike moments (GPZM) of order n and repetition m are defined as

$${\overline{Z}}_{nm}^{\alpha }={\int }_0^{2\pi }{\int }_0^1[{\overline{V}}_{nm}^{\alpha }(r,\theta ){]}^{\ast }f(r,\theta )rdrd\theta ,$$where * denotes complex conjugate. The orthogonality property of the pseudo-Zernike polynomials implies that

$$\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{1}{\int }}{\overline{V}}_{nm}^{\alpha }(r,\theta )[{\overline{V}}_{kl}^{\alpha }(r,\theta ){]}^{\ast }rdrd\theta ={\delta }_{nk}{\delta }_{ml}.$$

The image can be approximately reconstructed using the inverse transform formula (Teague, 1980)

$$f(r,\theta )\approx \sum _{n=0}^{\mathrm{\infty }}\sum _{|m|\le n}{\overline{Z}}_{nm}^{\alpha }{\overline{V}}_{nm}^{\alpha }(r,\theta ).$$

Bayesian estimator of classification accuracy

Let X_i, i = 1,2,…,s be the number of correctly predicted samples for the ith species (up to s species). These are the diagonal entries of the s × s classification matrix. For a statistical learning model, assume that it has constant probability π of correctly predicting the species identity of a sample (i.e. accuracy). Then, X_i is binomially distributed with number of trials equal to the number of i-th species in the sample (n_i) and success probability π. By Bayes’ Theorem, the posterior distribution of π, given X₁,X₂,…,X_s, is

$$f(\pi |X_1,X_2,\dots ,X_s)=P(X_1,X_2,\dots ,X_s|\pi )f(\pi )\times {({\int }_0^{1}P(X_1,X_2,\dots ,X_s|\pi )f(\pi )d\pi )}^{-1},$$where f(π) is the prior distribution of π. Using the conservative uniform prior f(π) = 1, 0 < π < 1, and assuming that $P(X_1,X_2,\dots ,X_s|\pi )=\prod _{i=1}^{s}P(X_i|\pi )$, it can be shown that f(π|X₁,X₂,…,X_s) has a beta distribution with shape parameters $\alpha =\sum _{i=1}^s{X}_i+1$, and $\beta =N-\sum _{i=1}^s{X}_i+1$, where N is the total sample size. Thus, the Bayesian posterior mean estimate of π is given by

$$\hat{\pi }=\frac{\sum _{i=1}^s{X}_i+1}{N+2}.$$

The lower and the upper limit of the 95% Bayesian credible interval of π are computed as the 2.5th and the 97.5th percentile of the beta distribution with shape parameters $\alpha =\sum _{i=1}^s{X}_i+1$ and $\beta =N-\sum _{i=1}^s{X}_i+1$, respectively.