Extreme learning machine for identifying soil-dwelling microorganisms cultivated on agar media
Sustained growth of yields is needed to fulfill the rise of consumption associated with increasing human population. For this purpose, scientists turn their attention towards the analysis of soil microorganisms that affect plant crop. The influence of pathogenic and other deleterious microorganisms could be very severe and may lead to enormous losses . In contrast, some species of soil microorganisms can effectively control plant pathogens and stimulate plant growth increasing amount and quality of crop . Understandably, the identification of microorganisms is critical in this area of research. Conventional methods are based on isolating the microorganism samples from soil, observing them under a microscope and lastly identifying examined organisms on the basis of their macro and micromorphology combined with phenotypic and molecular biological methods . The goal is to alleviate the latter by automatizing the task of soil microorganisms recognition based only on morphological traits. Noticeably such identification setting is only possible for some soil microorganism species and only on genera level due to the small differences in appearance of microorganisms within a genera . For this reason recognition conducted by microbiologists is a process that involves many different methods combined together that are time and cost consuming . Extensive studies have explored the application of Machine Learning (ML) for detecting pathogenic fungi based on leaf images, achieving remarkable accuracy levels . Convolutional Neural Network (CNN) Models such as AlexNet, GoogleNet, InceptionV3 and ResNets reach AUC scores close to 99% in detecting fungal pathogens through the images of leaves , . These studies underlie the potential of ML to accurately identify plant diseases and to classify the specific genera of microorganisms responsible for these state. Unfortunately, such approaches identify threats only after symptoms have been manifested. The resulting delay can lead to significant losses, especially with diseases like Verticillium wilt . It is particularly insidious as it spreads rapidly and can persist in the soil for extended periods, making it difficult to eradicate once it becomes established . Thus, early detection is crucial in managing such diseases before they cause widespread damage and agricultural losses. There are various time-consuming and costly methods for identifying microorganisms in soil samples, including phenotypic and molecular techniques . An alternative approach is to identify the genera of microorganisms with the aid of microscopic imaging. Noteworthy, while identifying soil fungi remains a niche study area, some insights might be gained from studies focusing on infections that are hazardous to human health. Zielinski et al. explored a similar challenge, demonstrating that the use of microscopic images can accelerate the identification process . Their research shows that combining CNNs with bag-of-words techniques can greatly reduce both the time and cost of fungal species identification by removing the necessity for extended biochemical tests. However, the limitation of Zielinski’s work is its reliance to the relatively small dataset of 180 images. Additionally, they all are captured under consistent lighting conditions and using a single type of microscope, with expert involvement potentially introducing biases. In contrast, the method presented in this paper not only addresses the need for rapid and cost-effective identification of soil microorganisms but also establishes a robust testing framework. The proposed model is validated across multiple datasets, each prepared under various conditions. This represents also the test scenario under which dataset is generated fully automatically, with the microscope capturing images sequentially and non-overlappingly to cover the sample area. Recent research has increasingly indicated the transformative potential of ML, particularly deep learning, in microbiology. In particular, studies by Treebupachatsakul and Khasim showcase the application of deep learning for recognizing and classifying microorganisms, with Khasim specifically emphasizing the efficacy of CNNs. Qu and Jiang provide comprehensive overviews of ML applications in microbiology, underscoring its utility in classification tasks and its potential to enhance the organization and application of microbiological knowledge. The surge in research utilizing CNN methods is notable due to their straightforward processing pipeline, which eliminates the need for labor-intensive feature crafting and often bypasses feature selection entirely. Despite this, traditional machine learning methods that involve feature engineering, such as Support Vector Machines (SVM), Random Forests (RF) and k-Nearest Neighbors (KNN), remain prevalent in microbial studies. For instance, Kotwal discusses their use in bacterial classification, demonstrating their continued relevance. Rani , in a comprehensive review covering 100 studies on machine learning applications in microbial recognition from 1995 to 2021, notes that only 12.1% of these papers focus on fungi. This statistic highlights the novel contribution of the current paper in addressing the identification of fungal microorganisms. An example of such focused research is the work by Liu et al. , which employs ML techniques to detect and count fungal microorganisms in microscopic images. Similarly, Tahir et al. utilized SVM for detecting fungal spores. Their approach involved extracting image patches, preprocessing them with Gaussian filters and using handcrafted features to achieve a notable accuracy of 88%. These studies collectively emphasize the significant and growing role of ML in advancing the field of microbiology, particularly in the area of fungal identification. The main research question explored in this study is whether a robust ML model can be developed for the accurate identification of specific soil microorganism genera across various datasets, including those generated through automated image capture. This capability is crucial for advancing systems that provide rapid and precise identification, thereby enabling timely interventions in agriculture to mitigate the devastating impacts of microorganism pathogens on crops or laboratory colonies. The datasets used in this research is obtained from Symbio-Bank - the collection of microorganisms collected by The National Institute of Horticultural Research in Skierniewice. A system developed here is designed to classify microscopic images of the following soil microorganisms - fungi of genera Fusarium , Trichoderma , Verticillium , Purpureolicillium and Chromista of genera Phytophthora . Only a single type of microorganisms is assumed to be present on the image upon applying rectification process to obtain monocultural isolate from the input pictures of the examined material. Microscopic images of soil-dwelling microorganisms are very complex and difficult to handle by standard image processing methods . Microorganisms exhibit a wide range of shapes and textures, including non-rigid and irregular forms, which complicates segmentation and analysis , . The respective organisms’ structures are not homogeneous as in the case of bacteria . Indeed, they are featured by various irregular structures like hyphae, phialides, micro and macroconidia, conidia cells, zoospores and many more. Watanabe conducted a brief study providing microbiologists a resource for rapid comparison of microorganism morphologies. Taxonomic differences among types of soil microorganisms include different shapes and morphology of microorganism/spore cells, different thickness of microorganism/spore cell walls, differences in color, microorganism/spore cell structure. Consequently, the goal of constructing a classifier using machine vision methods constitutes a non-trivial task. The research leverages five distinct datasets specifically curated for this study. The initial dataset, comprising 128 images, was meticulously prepared by an expert microbiologist and served as a basis for fine-tuning the image processing pipeline. This foundational set was expanded into a second dataset with 303 images. Subsequently, three additional datasets were created using automated image acquisition techniques, collectively contributing to a total of 2,866 images. During the literature review, it was observed that this research employs the largest database in publications dedicated to soil microorganism identification. Typically, manual image acquisition is a labor-intensive process, leading most studies to work with datasets containing fewer than 500 images. Furthermore, these images are often captured under a single setup and methodology, which may introduce challenges in generalizing the technology to new environments. To address these limitations, the model was initially optimized using the first dataset, then trained and validated across four additional datasets. To ensure the findings are both robust and interpretable, the study integrates Explainable AI (XAI) techniques. This approach not only evaluates the model’s performance but also elucidates the rationale behind the classification of individual instances and overall test sets, as elaborated in the Results section, which also includes standard performance metrics. A key aspect of the study is also the comparison between traditional handcrafted feature approaches and CNNs. On the other hand, study aims to advance the research towards practical application stages, where prediction time is critical, especially during large-scale tests. For this purpose, the study employs the Extreme Learning Machine (ELM), a type of Single Layer Feedforward Neural Network (SLFNN), known for its rapid training and exceptionally fast prediction times. ELMs are particularly suited for practical applications where quick predictions are essential . The presented system can distinguish only the tested types of soil microorganisms based on spores obtained from cultures of microorganisms on agar media and visualized on microscopic images. Within the genera of soil microorganisms, there is variability at the level of species and strains , . At this stage of the development of the presented system, the authors focused on the genera of microorganisms. The authors realize this method will fail flawlessly due to phenotypic variability within species and genera. This paper accurately recognizes the selected soil-dwelling using the ML model that, in this stage, can be directly used in the contamination early-detection of laboratory cultures by the analyzed genera of microorganisms. Detection of microorganisms in soil directly on the soil photo is currently not possible, and up to now, any system capable of performing such a task has not been proposed; this is the scope for future development merging methods proposed in this paper of phenotypic recognition with molecular as the ultimate system. For isolates belonging to Fusarium , Verticillium , Phythophthora , Trichoderma and Purpureocillium , which do not grow on agar media, it will not be possible to use the presented method. The demonstrated approach is the first stage of work aimed at differentiating among themselves, which is critical for horticultural practice in the five types of soil microorganisms.   In this section, the results from various evaluations are comprehensively presented, focusing on the detailed analysis of the features generated across different datasets. The primary classifier employed in this research is the ELM. To provide a comprehensive evaluation, comparisons with different ML algorithms are also presented in terms of classification metrics and time of execution. Feature analysis This subsection explores the distinctions and commonalities in the features across the various datasets. Each dataset comprises a set of 74 features, as detailed in the Methods section. To initiate this analysis, boxplots are utilized to visually represent the value distributions of specific features across the different datasets (see Fig. ). This method highlights the statistical differences among datasets by selecting representative features from each category for visualization. The boxplots effectively illustrate how the statistical properties of features vary between datasets, underscoring the diversity in the characteristics of images from different strains, microscopes and lightning conditions. Multivariate Analysis of Variance (MANOVA) was employed to assess the differences between datasets across a reduced feature space (10 dimensions) derived from Principal Component Analysis (PCA) . The MANOVA revealed significant findings across both the intercept and dataset comparisons. In terms of the intercept, the analysis indicated that the combined effect of all PCA components was significantly different from zero (Wilks’ lambda = 0.3304, F(10, 2818) = 571.1673, p< 0.0001). This suggests a robust multivariate separation in the reduced feature space. Furthermore, comparing datasets showed highly significant differences (Wilks’ lambda = 0.0119, F(40, 10687) = 592.4963, p < 0.0001), indicating distinct multivariate distributions among the datasets. Pillai’s Trace (2.5217) shows a significant portion of variance in the dependent variables is explained by the grouping variable (Dataset label). Hotelling-Lawley Trace (9.6309) suggests substantial multivariate differences between the datasets. Roy’s Greatest Root (5.1398) implies a strong difference in the first canonical dimension underscoring the statistically meaningful differences of feature values across the datasets. Visual confirmation through t-Distributed Stochastic Neighbor Embedding (t-SNE) analysis in Fig. further supported these findings, illustrating clear separation between datasets in the reduced-dimensional space highlight the different nature of the datasets employed in this research yielding the complex task for ML model. Performance of the model The performance of the ELM is evaluated both on individual datasets and on combined datasets (marked as All) to assess its robustness and adaptability. The ELM performance is tested on individual datasets, varying the number of neurons in the hidden layer from 100 to 2000. The activation function used is Mish, an unconventional choice for ELM but one that recent studies suggest it can be successfully applied to ELM, borrowing from its use in deep learning techniques . The implementation utilized is TfELM , which accelerates computations by leveraging the TensorFlow library to run on a GPU. This setup enables efficient processing and enhances the ELM’s computational performance. The evaluation follows a protocol of 50 iterations of 10-Fold Cross-Validation to ensure reliable and stable results for a given dataset. The ultimate performance metrics, representing the best outcomes for the specified number of neurons, are summarized in Table . The ROC/AUC scores reflect the one-vs-rest approach used for multi-class classification. The results indicate exceptional performance, as evidenced by high F1 scores, precision, recall and near-perfect ROC AUC values across all datasets, affirming the model’s robustness and accuracy in classifying various genera of soil microorganisms by the accurate feature generation from various datasets. Comparison to the other models The performance of ELM was benchmarked against several well-established classifiers to evaluate its efficiency and accuracy. The chosen classifiers for this comparison included Multi-Layer Perceptron (MLP), Gradient Boosting, XGBoost, LightGBM, CatBoost and Random Forest , all implemented using the scikit-learn library. ELM demonstrated promising results, particularly comparable to CatBoost, highlighting its capability for precise predictions. Table represents the comparison tests across different datasets and methods. For reference, experiments were also conducted using Convolutional Neural Networks (CNNs) - ResNet50 and ResNet152V2, fine-tuned on pre-trained ImageNet models. Input images from each dataset were resized to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$224  224  3$$ and split into training, validation, and testing sets in a 70:20:10 ratio. The models were adapted with a custom top layer consisting of two fully connected layers with 512 neurons, separated by a Dropout layer with a rate of 0.3 to nullify 30% of the weights, enhancing robustness and mitigating overfitting. The Adam optimizer with a learning rate of 0.004 and the Sparse Categorical Crossentropy loss function were used. Each model was trained for a maximum of 1000 epochs, with Early Stopping implemented to halt training if the validation loss did not decrease further. The results are summarized in Table . Time execution The computational efficiency of various classifiers was assessed to evaluate their suitability for the system. The evaluation involved measuring both the average fit time (training duration) and the average prediction time during a single fold of 10-Fold cross-validation. For this purpose, the largest dataset was used, a combination of Dataset1-5, comprising 2,866 images. The results, summarized in Table , clearly illustrate the exceptional performance of the ELM in terms of both fit and prediction times. Compared to other classifiers, ELM stands out, particularly when contrasted with CatBoost that is proved to provide comparable accuracy performance with ELM. The ELM demonstrated a significant advantage, being six times faster in the fit phase and twice as fast in the prediction phase, underscoring its potential for practical industrial applications. Given its accuracy and speed, ELM offers substantial benefits over other classifiers. Noteworthy, CNN time of the prediction has been measured as 0.11 s, twelve times longer than ELM. Explainable artificial intelligence This subsection explores the transparency and interpretability of the model’s outputs by employing XAI techniques. To elucidate the impact of the handcrafted features, Shapley Additive Explanations (SHAP) are utilized. SHAP values provide insights into the contribution of each feature to the model’s predictions for any machine learning method . Figure presents a summary plot for Dataset1, where the ELM with 1000 neurons in the hidden layer and the Mish activation function is analyzed. This plot allows us to discern whether specific features have a positive or negative influence on the model’s predictions, thereby enhancing the interpretability of the classification decisions. The beeswarm plot provides a clear visualization of the most influential features in the overall classification process using the test set. As observed, a variety of feature types contribute significantly, indicating that these features serve as robust descriptors of image characteristics. Additionally, individual prediction explanations can be examined using the waterfall plot, which illustrates how specific traits positively or negatively influence the model’s prediction for a given sample (see Fig. ). This dual approach-examining both the overall feature importance and the detailed contribution for individual samples-enhances the understanding of the model’s decision-making process. Notably, the analysis of both beeswarm and waterfall SHAP plots reveals that the selected feature sets, including statistical measures over the I1I2I3 and XYZ color spaces, as well as GLCM and GLRLM features, each contribute significantly to the final decisions made by the classifier. Another approach involved using SHAP to evaluate the overall influence of features on the target variable. The Dataset1 was cleaned, normalized, and split into training and testing subsets (80/20). An ELM model with 1,000 neurons and a mish activation function was trained. SHAP explainability was applied to the testing data to calculate SHAP values and assess feature importance. The absolute SHAP values were summed across all samples, allowing the ranking of features based on their contribution to the model’s predictions. The most influential features were identified and are presented in Table .   This subsection explores the distinctions and commonalities in the features across the various datasets. Each dataset comprises a set of 74 features, as detailed in the Methods section. To initiate this analysis, boxplots are utilized to visually represent the value distributions of specific features across the different datasets (see Fig. ). This method highlights the statistical differences among datasets by selecting representative features from each category for visualization. The boxplots effectively illustrate how the statistical properties of features vary between datasets, underscoring the diversity in the characteristics of images from different strains, microscopes and lightning conditions. Multivariate Analysis of Variance (MANOVA) was employed to assess the differences between datasets across a reduced feature space (10 dimensions) derived from Principal Component Analysis (PCA) . The MANOVA revealed significant findings across both the intercept and dataset comparisons. In terms of the intercept, the analysis indicated that the combined effect of all PCA components was significantly different from zero (Wilks’ lambda = 0.3304, F(10, 2818) = 571.1673, p< 0.0001). This suggests a robust multivariate separation in the reduced feature space. Furthermore, comparing datasets showed highly significant differences (Wilks’ lambda = 0.0119, F(40, 10687) = 592.4963, p < 0.0001), indicating distinct multivariate distributions among the datasets. Pillai’s Trace (2.5217) shows a significant portion of variance in the dependent variables is explained by the grouping variable (Dataset label). Hotelling-Lawley Trace (9.6309) suggests substantial multivariate differences between the datasets. Roy’s Greatest Root (5.1398) implies a strong difference in the first canonical dimension underscoring the statistically meaningful differences of feature values across the datasets. Visual confirmation through t-Distributed Stochastic Neighbor Embedding (t-SNE) analysis in Fig. further supported these findings, illustrating clear separation between datasets in the reduced-dimensional space highlight the different nature of the datasets employed in this research yielding the complex task for ML model.   The performance of the ELM is evaluated both on individual datasets and on combined datasets (marked as All) to assess its robustness and adaptability. The ELM performance is tested on individual datasets, varying the number of neurons in the hidden layer from 100 to 2000. The activation function used is Mish, an unconventional choice for ELM but one that recent studies suggest it can be successfully applied to ELM, borrowing from its use in deep learning techniques . The implementation utilized is TfELM , which accelerates computations by leveraging the TensorFlow library to run on a GPU. This setup enables efficient processing and enhances the ELM’s computational performance. The evaluation follows a protocol of 50 iterations of 10-Fold Cross-Validation to ensure reliable and stable results for a given dataset. The ultimate performance metrics, representing the best outcomes for the specified number of neurons, are summarized in Table . The ROC/AUC scores reflect the one-vs-rest approach used for multi-class classification. The results indicate exceptional performance, as evidenced by high F1 scores, precision, recall and near-perfect ROC AUC values across all datasets, affirming the model’s robustness and accuracy in classifying various genera of soil microorganisms by the accurate feature generation from various datasets.   The performance of ELM was benchmarked against several well-established classifiers to evaluate its efficiency and accuracy. The chosen classifiers for this comparison included Multi-Layer Perceptron (MLP), Gradient Boosting, XGBoost, LightGBM, CatBoost and Random Forest , all implemented using the scikit-learn library. ELM demonstrated promising results, particularly comparable to CatBoost, highlighting its capability for precise predictions. Table represents the comparison tests across different datasets and methods. For reference, experiments were also conducted using Convolutional Neural Networks (CNNs) - ResNet50 and ResNet152V2, fine-tuned on pre-trained ImageNet models. Input images from each dataset were resized to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$224  224  3$$ and split into training, validation, and testing sets in a 70:20:10 ratio. The models were adapted with a custom top layer consisting of two fully connected layers with 512 neurons, separated by a Dropout layer with a rate of 0.3 to nullify 30% of the weights, enhancing robustness and mitigating overfitting. The Adam optimizer with a learning rate of 0.004 and the Sparse Categorical Crossentropy loss function were used. Each model was trained for a maximum of 1000 epochs, with Early Stopping implemented to halt training if the validation loss did not decrease further. The results are summarized in Table .   The computational efficiency of various classifiers was assessed to evaluate their suitability for the system. The evaluation involved measuring both the average fit time (training duration) and the average prediction time during a single fold of 10-Fold cross-validation. For this purpose, the largest dataset was used, a combination of Dataset1-5, comprising 2,866 images. The results, summarized in Table , clearly illustrate the exceptional performance of the ELM in terms of both fit and prediction times. Compared to other classifiers, ELM stands out, particularly when contrasted with CatBoost that is proved to provide comparable accuracy performance with ELM. The ELM demonstrated a significant advantage, being six times faster in the fit phase and twice as fast in the prediction phase, underscoring its potential for practical industrial applications. Given its accuracy and speed, ELM offers substantial benefits over other classifiers. Noteworthy, CNN time of the prediction has been measured as 0.11 s, twelve times longer than ELM.   This subsection explores the transparency and interpretability of the model’s outputs by employing XAI techniques. To elucidate the impact of the handcrafted features, Shapley Additive Explanations (SHAP) are utilized. SHAP values provide insights into the contribution of each feature to the model’s predictions for any machine learning method . Figure presents a summary plot for Dataset1, where the ELM with 1000 neurons in the hidden layer and the Mish activation function is analyzed. This plot allows us to discern whether specific features have a positive or negative influence on the model’s predictions, thereby enhancing the interpretability of the classification decisions. The beeswarm plot provides a clear visualization of the most influential features in the overall classification process using the test set. As observed, a variety of feature types contribute significantly, indicating that these features serve as robust descriptors of image characteristics. Additionally, individual prediction explanations can be examined using the waterfall plot, which illustrates how specific traits positively or negatively influence the model’s prediction for a given sample (see Fig. ). This dual approach-examining both the overall feature importance and the detailed contribution for individual samples-enhances the understanding of the model’s decision-making process. Notably, the analysis of both beeswarm and waterfall SHAP plots reveals that the selected feature sets, including statistical measures over the I1I2I3 and XYZ color spaces, as well as GLCM and GLRLM features, each contribute significantly to the final decisions made by the classifier. Another approach involved using SHAP to evaluate the overall influence of features on the target variable. The Dataset1 was cleaned, normalized, and split into training and testing subsets (80/20). An ELM model with 1,000 neurons and a mish activation function was trained. SHAP explainability was applied to the testing data to calculate SHAP values and assess feature importance. The absolute SHAP values were summed across all samples, allowing the ranking of features based on their contribution to the model’s predictions. The most influential features were identified and are presented in Table .   The central objective of this study was to develop a robust ML pipeline capable of accurately classifying soil microorganisms based on their microscopic images. This task, particularly focusing on specific genera of fungi and Chromista , presents significant challenges due to the inherent complexity and variability of these organisms’ structures. The results from the extensive experiments underscore the proposed approach effectiveness highlighting substantial advancements in automated microorganism identification. Across all datasets, ELM consistently achieved high accuracy, precision, recall and ROC/AUC values, demonstrating its robustness and adaptability to varying imaging conditions and microorganism strains. This consistency is notable whether the model was tested on individual datasets or a combined dataset, underscoring its potential for practical application in diverse environments. When benchmarked against other well-established machine learning models, including MLP, Gradient Boosting, XGBoost, LightGBM, CatBoost and RF, ELM outperformed these methods in terms of both speed and accuracy. Its superior balance of computational efficiency and predictive power makes it particularly suitable for large-scale data processing tasks. The comparison of classical feature-based models and CNNs highlights the performance of ELM in terms of execution time and predictive accuracy. While the analyzed CNN models generally achieved higher F1 scores and AUC values across most datasets, their prediction time (0.11s) was significantly slower compared to ELM (0.009s). This makes ELM particularly advantageous in scenarios where quick predictions are critical. Despite its lower performance metrics in some cases, ELM demonstrated competitive results, achieving an F1 score of up to 0.984 in specific datasets, showcasing its balance between computational efficiency and predictive power. The comprehensive analysis of the datasets using MANOVA revealed significant multivariate differences among them. This finding indicates that each dataset captures unique aspects of the microorganism images, emphasizing the importance of using a diverse set of data to train and validate the model. The variability among datasets is crucial for ensuring that the model is generalizable and robust across different imaging conditions and microorganism characteristics. Another contribution of this study is the use of SHAP values to understand the impact of individual features on the model’s predictions. This analysis provides crucial insights into how the model interprets different morphological traits and makes classification decisions. The ability to interpret these decisions is vital, especially in scientific and industrial applications where understanding the model’s reasoning can enhance trust and acceptance of its outputs. The most influential features, as identified by SHAP, span across all categories of the generated features described in the Feature Calculation section. This highlights the effectiveness of the approach, with texture features, in particular, emerging as the most descriptive for the ELM model. The proposed approach also incorporated datasets obtained through automated image acquisition. This automation not only expanded the dataset size significantly but also introduced greater variability in imaging conditions. The model’s successful performance on these automated datasets demonstrates its potential for deployment in real-world diagnostic systems, where automated and rapid identification processes are essential. Noteworthy, in particular ELM performed better for datasets automatically acquired than the manual one taken by the expert. Comparing research with existing literature, it becomes clear that this study sets a new benchmark in the field of soil microorganism identification. The use of a substantially larger and more varied dataset, coupled with automated image acquisition and a robust machine learning model, juxtapose this work from earlier efforts that often relied on smaller, manually curated datasets and more uniform imaging conditions. Future research could focus on scaling this approach for broader use, including integration into automated imaging systems for early detection of pathogenic microorganisms. Additionally, while this study focused on specific genera of fungi and Chromista , the methodology could be expanded to include a broader range of soil microorganisms, such as bacteria, by addressing their unique morphological characteristics. A key challenge in such an expansion would involve crafting features tailored to the new genera. Adding new genera to the analyzed dataset would require evaluating the existing features to ensure they effectively capture the relevant morphological differences. While the current feature set emphasizes color and texture, new datasets that include genera with distinct structural traits could benefit from the integration of additional geometric or shape descriptors, such as circularity, elongation, or convexity. Moreover, incorporating multi-scale texture features could improve the detection of more complex structures, complementing the existing color and texture-based metrics. If new genera are imaged using alternative methods, such as fluorescence microscopy, the system could be adapted by recalibrating color space features or incorporating intensity-specific metrics tailored to these imaging conditions. Combining the machine vision-based approach with traditional microbiological methods could yield a hybrid system that leverages the strengths of both. For instance, rapid initial screening using proposed model could be followed by targeted molecular or phenotypic tests for confirmatory diagnosis, thus enhancing both speed and accuracy. Another potential application is the quick detection of cross-contamination during the cultivation process. This is particularly critical when working with soil-dwelling microorganisms, where multiple genera coexist and interact, potentially affecting experimental outcomes and product quality. Contamination detection is especially challenging with certain genera, such as Trichoderma , Fusarium , Verticillium , and Purpureocillium . A hybrid system combining rapid machine vision analysis with confirmatory microbiological methods could be invaluable for detecting and mitigating such contamination, ensuring the integrity of experimental cultures. To ensure the model’s adaptability to new and diverse environmental conditions, future studies should explore its performance on datasets obtained from different geographic regions, soil types and agricultural practices. This would further test the model’s robustness and ensure its applicability in a wide range of real-world scenarios. In conclusion, this study successfully addresses the complex task of classifying soil microorganisms based on microscopic images. By developing and validating a robust ML model and integrating diverse datasets including those obtained through automated image acquisition, the advancement in the field is significant. The insights gained from feature analysis and the application of XAI techniques underscore the model’s utility and transparency. This work not only advances automated microorganism identification but also lays the groundwork for future innovations in microbial diagnostics and soil health monitoring. The proposed approach holds great promise for transforming agricultural practices through early detection and management of soil-borne pathogens, contributing significantly to sustainable crop production and food security.   In this study, a comprehensive pipeline is proposed, integrating image preprocessing, segmentation, feature generation and classification using extracted features. Each element of this pipeline is meticulously elaborated upon in this section. Additionally, the five distinct datasets employed in this research are thoroughly explained, along with the detailed processes of image preprocessing, segmentation, feature extraction and the specifics of the ELM method. Datasets Microorganisms are grown in Potato Dextrose agar, specifically Merck 1.10130.0500 at 26 degrees Celsius. In the case of Fusarium , Verticillium and Phytophthora the microorganisms are incubated for 7–10 days, while for Trichoderma strain it takes about 4-5 days until conidial spores are produced. In total five different datasets were formed that can be obtained from the corresponding author upon the reasonable request. Dataset1 The initial experiments utilized a dataset comprising 128 images of microorganisms from the genera Verticillium , Trichoderma , Phytophthora and Fusarium , referred to as Dataset1 throughout this manuscript with exemplary images presented in Fig. . This foundational dataset is based on the work of Kruk et al. , which demonstrated the potential for classifying soil-dwelling microorganisms using ML techniques. The primary objective with Dataset1 was to develop a robust image preprocessing and segmentation pipeline that could accurately differentiate the microorganism regions from the background as this is recognized as a crucial step for ultimate model’s performance. Despite being well-suited for preliminary investigations into image segmentation and feature extraction, Dataset1 was relatively small, which limited its ability to ensure the robustness of the solution across diverse scenarios. To overcome these limitations and to provide a comprehensive evaluation of the methodology’s performance, the additional, substantially larger datasets are provided. These new datasets were crucial in validating the generalizability and effectiveness of the proposed image preprocessing, masking and classification methods at scale. Dataset1 was carefully curated by professional microbiologists who were responsible for both sample preparation and image acquisition through microscopy. The primary goal was to generate high-quality images that clearly showcased the morphological characteristics essential for accurate microorganism identification. To achieve this, each image was meticulously focused to ensure that the unique features of each genus were prominently visible. The process began with the cultivation of the microorganisms, followed by the careful transfer biological material onto microscope slides. To facilitate optimal imaging, a drop of distilled water was added to dilute the biological material and a cover slip was placed on top (see Fig. ). These prepared slides were then subjected to microscopic examination, allowing the microbiologists to capture detailed images at different magnification levels that would serve as the basis for the ML analysis. Given that the model in question solely relies on image data for input without the incorporation of additional multimodal features, the quality and clarity of these images were paramount. However, this approach introduces a degree of bias, as image acquisition process rely heavily on the expertise of the microbiologists. This necessitates the involvement of experts throughout the pipeline to ensure that the images accurately represent the microorganisms’ traits for effective model training and evaluation. Dataset2 Dataset2 builds upon the foundation of Dataset1 by expanding its scope and content. This expanded dataset includes additional images of the original genera ( Trichoderma and Fusarium ) and introduces images of a new genus, Purpureocillium (see Fig. ). Prepared using the same procedures as Dataset1, Dataset2 effectively consolidates and enhances the original set, resulting in a comprehensive collection of 304 images. Dataset3 and Dataset4 Upon evaluating Dataset1 the accuracy of the generated segmentation masks is assessed for the suitability of the calculated features and the overall performance of the applied ML methods. This comprehensive analysis highlighted the need for a larger, more automated dataset to minimize bias and enhance scalability. Thus, the third dataset, referred to in this paper as Dataset3 was prepared, aiming to automate the entire image acquisition process and build a significantly larger dataset without the direct involvement of experts. The critical point in creating Dataset3 was the elimination of expert bias inherent in manually capturing images, which, while valuable for understanding morphological traits, does not reflect real-world scenarios where a fully automated system is the target. The approach involved using an automated microscope equipped with a movable mechanical stage controlled by built-in servo motors and software drivers. This setup allowed the microscope to systematically capture images across a predefined area without overlap. During initial trials, the challenges with out-of-focus images due to the varying depths (Z-axis) of microorganisms at high magnification levels and slight imperfections on the slide surfaces were encountered. Although, the autofocus feature helped address this issue, it occasionally failed to achieve perfect focus, leading to some variability in image quality. Despite these challenges, the automated system substantially reduced the reliance on expert involvement and increased the efficiency of image acquisition. To ensure comparability, the samples for Dataset3 were cultivated similarly to those in Dataset1. Two methods of material transfer onto microscope slides were examined. The first method, consistent with Dataset1, involved placing the biological material on a slide, adding a drop of water and covering it with a slip. This approach created Dataset3. The second method involved smearing the material directly onto the slide without a cover slip, which formed a new dataset, referred to as Dataset4 (see Fig. b). Both Dataset3 and Dataset4 were collected without requiring biological expertise for each image. Instead, a trained man handled the sample preparation and setup of the microscope, which included positioning the sample on the mechanical stage, using software to target areas of interest and initiating the automated imaging process with autofocus. This approach simulates a real-world scenario where minimal expert intervention is needed. Sample images from both datasets are presented in Fig. . Looking forward, there is potential for further automation, such as using robots for material transfer and advanced software for precise targeting, which would enhance the efficiency and reliability of the image acquisition process. These advancements are crucial for the ongoing development and scalability of the system. Dataset5 The fifth dataset, referred to as Dataset5, was generated following a procedure similar to that of Dataset4 (Fig. ). For Dataset4, images were captured in a non-overlapping, systematic manner to cover the entire sample area, as illustrated in Fig. b. Each image was taken using the microscope’s autofocus feature, which aimed to optimize focus for each shot. However, while autofocus provided a reasonable solution, it occasionally failed to achieve perfect focus across all images. To address the limitations of autofocus, Dataset5 incorporated a depth imaging technique. This method involved capturing a series of images at varying focal distances (the distance between the sample and the camera lens) for each section of the sample. These multiple images, each focused at a different depth, were then automatically stacked to produce a single composite image. This stacking process significantly enhances the overall image quality by ensuring that all areas of the sample are in focus, thus mitigating the problem of focal inconsistencies observed with single-focus images. Although this depth imaging approach enhances image clarity significantly, it requires a slower data acquisition process due to the need to capture and stack multiple images at different focal depths for each area of the sample. For Dataset5, the additional challenges were encountered when attempting to use the depth imaging technique in conjunction with the traditional method of adding a drop of water between the microscope slide and the coverslip that was applied for Dataset3. This method, commonly used to stabilize the sample for conventional imaging, proved unsuitable for depth imaging. The water layer caused instability, leading to slight movements of fungal elements within the sample. Consequently, these movements resulted in blurred or misaligned images when stacking the focal series, thereby undermining the intended clarity of the depth imaging process. Overview of the datasets Initially, Dataset1 was prepared for preliminary research, focusing on data segmentation and feature calculation. This dataset included 128 images from four genera: Verticillium , Trichoderma , Phytophthora and Fusarium . To test and to enhance the robustness of the solution, this dataset was expanded by creating Dataset2. This expanded dataset added more images for Trichoderma and Fusarium and incorporated a new genus, Purpureocillium , increasing the total to 303 images. Subsequently, with the implementation of an automated image acquisition system, three additional datasets were developed: Dataset3, Dataset4 and Dataset5. These datasets were captured using an automated microscope, allowing to gather a significantly larger and more diverse collection of images. Unfortunately, Phytophthora could not be included in the automated datasets due to deterioration of all samples held prior to the establishment of the automated imaging system. Collectively, these efforts resulted in the creation of the largest database of soil microorganism images captured through microscopy for the ML purposes, setting a new standard in the literature. Datasets with all details are provided in Table . The findings presented in this research are crucial due to the successful replication of the entire machine vision algorithm and ML methods across a range of datasets. These datasets vary significantly in size and were acquired using both automated and manual methods, as well as through two distinct types of microscopes. Such replication demonstrates that the system retains its performance consistently across different conditions, indicating a strong resistance to overfitting. For statistical proofs of this differences please refer to section Feature analysis . This robustness is evident from the ability to maintain high accuracy and reliability in image analysis and classification, regardless of how the data was collected. The diverse operational settings in which the datasets were obtained, including automated and non-automated means, further validate the adaptability of the system. This versatility ensures that the system may be successfully deployed in a variety of real-world circumstances where consistency and reliability are critical. By proving effective across a broad spectrum of datasets and acquisition techniques, the system showcases its potential for scalability and generalization in the field of microbial image analysis. Image preprocessing All input images are automatically preprocessed in order to improve the segmentation phase which in turn has an essential impact on final classification accuracy. Indeed, the primary challenge associated with microscopic images lies in managing noise interference and low contrast impact. Another significant consideration is the presence of scale labels positioned at the bottom-left corner of each image for Dataset1 and Dataset2 only. It is important to note that the colors of these labels may vary which potentially brings distortion to the generated results. The latter occurs due to the automatic color selection of the labels by a microscopic software. Image masks that differentiate image to label image foreground and background are prepared by thresholding, whose values are empirically chosen. Consequently, the labels are concealed by replacing their area with the average color value of the complete image. The following stage involves converting the given image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}$$ in RGB color space to grayscale forming image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{gray}$$ , which can be subjected to further processing. To mitigate image noise, this study applies Gaussian and Wiener filters . Finally, to improve the clarity of image features and to aid in distinguishing microorganisms from the background, the Contrast Limited Adaptive Histogram Equalization (CLAHE) algorithm is employed. Image segmentation Majority of the microscopic images in the datasets exhibit an imbalanced ratio between the area of the background and the area of the region of interest (ROI) that contains the soil microorganisms. Most of the area on microscopic pictures is covered by background information irrelevant to the classification task e.g. test-tube, medium for the cultivation of bacteria or minor impurities. Therefore, background has to be separated from the ROI. To address this issue, a binary image mask is created, assigning a value of 1 to each pixel that belongs to the ROI and a value of 0 to each pixel that belongs to the background. One of the most commonly used methods to perform image segmentation is an Otsu method that calculates one threshold value globally for the whole image upon minimizing weighted sum of variances within ROI and background pixel values and maximizing weighted sum of variances between them. According to , the Otsu method is known to produce weak results when the area of the ROI is imbalanced in comparison to the background area. Due to the above limitation a method called Adaptive Image Thresholding is used. The mask obtained using adaptive image threshold needs further image processing using morphological operations. Open operation is used to remove small objects in image and close operation is applied to join together objects that are close to each other. Both operations are applied to remove segmentation imperfections caused by the noise occurrence within original image (the so-called salt and pepper noise removal—see ). Despite incorporating morphological operations, the thresholding results remain unsatisfactory. Classification accuracy of the ELM classifier applied on features computed for the generated mask renders 82.7% for Dataset1. Unsatisfactory results are caused by the fact that input images are very complex as having different structures and textures. The segmentation accuracy is enhanced up to 87.4%, thus improving the performance of the classifier upon applying k -means algorithm . Mask M is multiplied using Hadamard product with [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ forming combined image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ . On that image k -means method with [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$k=3$$ is run. As a result the original image is separated into three classes: i.e. a background and ROI area subsequently divided to two disjoint subclasses. Pixels that are classified to the brightest centroid represent a new, smaller area of ROI. The pixels with the highest contrast are assumed to have the most relevant information to distinguish microorganisms. In contrast, ROI pixels assigned to the centroid with lower values are classified in fact to the background as it may contain data disrupting the final results. In this manner a new mask M is generated. Finally, on a given mask, open and close morphological operations are applied. All crucial image segmentation stjpg contributing to the above algorithm are illustrated in Fig. and the results from various datasets in Fig. . Features calculation Image features are computed as descriptors to measure statistical attributes of an image. The primary objective is to devise features that can effectively discriminate the different genera of microorganisms, enabling the formation of a robust classifier. Thus, color and texture-based handcrafted features are calculated in this study that are explained in the forthcoming subsections. Color features Color features are computed on combined images [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ generated by Hadamard product of original image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ with its mask M introduced in the previous section. Conversion of an image from RGB color values to a different color space prior to feature calculation can enhance classification accuracy . Indeed, experiments show such increase once conversion from RGB to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ color space is applied. First set of color features contains 8 statistical features, each calculated separately for three image color channels of color space [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ to which [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ was previously converted. Every image pixel of coordinates x and y in RGB image is a vector containing 3 values for the respective RGB color space i.e. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_{RGB}[x,y] = [R[x,y],G[x,y],B[x,y]]$$ . Conversion from RGB to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ is made upon applying the following operation: [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ I_{I_1I_2I_3}[x,y] = I_{RGB}[x,y] 1/3 & 1/3 & 1/3\\ 1/2 & 0 & -1/2\\ -1/4 & 1/2 & -1/4 . $$ Histograms for images [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ (denoted by H and h , respectively, generated upon conversion to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ color space) are calculated for every color space channel. Image histogram is a vector h , where h [ i ] is a number of occurrences of image pixels with the corresponding intensity value ranging within [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0 i 255$$ . Having both histograms the following features are calculated: mean and standard deviation of H; mean, standard deviation, kurtosis, skewness, entropy and energy of h . The second set of color features is based on image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ converted into XYZ color space. In contrast to the previous case, the calculated measures rely directly on image pixels, not histogram distribution. The latter permits to extract additional information on pixel spatial relationships which in turn is unavailable from image histogram . The XYZ color space is selected after comparison of classification accuracy for the different color spaces. For every X , Y and Z color channel four statistical measures are calculated separately i.e. mean, standard deviation, kurtosis and skewness yielding in total 12 features. Texture features Image texture defines spatial relations of chrominance or luminance of the image pixels. Grey-level co-occurrence matrix (GLCM) GLCM introduced by Haralick is a method to extract spatial relations of pixels for a given angle direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and maximal distance d . To apply the GLCM algorithm, the first step involves quantizing the image into q levels, resulting in the creation of an image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q$$ . Next, four co-occurrence matrices [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }_{q q}$$ are defined, each computed for a different angle [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ (i.e., 0, 45, 90 degrees and 135 degrees). Here all coefficients of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }_{q q}$$ are initialized to zero. Examination of pixels relations is performed over their local neighborhood using window W of size [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$m n$$ pixels. In sequel, W is iteratively moved throughout the image I assuming that coordinates of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q$$ are positioned in the top left corner. (see Fig. ). At the first iteration of algorithm matrix W contains pixels of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q[y,x]$$ with indices satisfying [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<y m$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<x n$$ . Next W is moved horizontally so that it does not overlap with the previous window - e.g. indices of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q[y,x]$$ for the second window range over [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<y m$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$n<x  2n$$ . If such left-to-right shift of W is not anymore possible such window is shifted vertically to underneath-left boundary position (e.g. after the first horizontal shifting movement the indices run over [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$m<y 2m$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<x n$$ as illustrated in Fig. (see window III). GLCM algorithm terminates once none of the above stjpg is further possible (i.e. W is too close to the right bottom corner of the I ). Every obtained W is the subject of further calculations. Co-occurrences of pixel values [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i=0, , q-1$$ are noted for a direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = 0^{ }, 45^{ }, 90^{ }, 135^{ }$$ and maximal distance d . Let [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{11}, ,w_{mn}$$ be coefficients of a given matrix W . Iterating through all of the coefficients one analyzes spatial correlations in their neighborhood. Currently analyzed pixel (coefficient of W in certain iteration) is named pixel of interest p . Definition of directions for maximal distance [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$d=2$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$p = w_{42}$$ is presented in Fig. . Assume p is [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{42}=0$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ =0^{ }$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$d=2$$ . Rows in [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ store the value of p and columns values of co-occurrence pixels being in relation to p . Since [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{42}=0$$ , the only coefficients in the first row of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ can be incremented. Then, based on Fig. , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{43}$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{44}$$ should be now analyzed. Coefficient [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{43}=0$$ implies that value in [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ in the first column ( [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}[1,1]$$ ) is incremented. Analogously, for [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{43}=3$$ [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}[1,4]$$ is increased. Both stjpg are presented in Fig. . In the same manner calculations for every W are performed analyzing all of the [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{ij}$$ for a direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ . Algorithm’s result is four matrices: [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{45^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{90^{ }}$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{135^{ }}$$ whose coefficients are then converted to the probabilities by dividing each of its values by mn . In the next step, the following statistical measures called GLCM features based on coefficients in matrices [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }$$ can be computed: Contrast, Correlation, Energy, Homogeneity, Autocorrelation, Cluster Prominence, Inverse Difference and Dissimilarity . Noticeably each of these features is calculable separately for each [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = 0^{ }, 45^{ }, 90^{ }, 135^{ }$$ . Having computed [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{0^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{45^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{90^{ }}$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{135^{ }}$$ their mean [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy = 0.25(Energy^{0^{ }}+Energy^{45^{ }}+Energy^{90^{ }}+Energy^{135^{ }})$$ is next derived. Originally GLCM based features are calculated on the grayscale image i.e. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_t^{gray}$$ . A further improvement is obtained here once GLCM features are computed separately for each of the RGB image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_t^{RGB}$$ channel yielding in total 24 features. In contrast to the work by Duan et al. , which proposes an extension that utilizes images in a different color space than RGB for GLCM feature computation, the experiments incorporating different color spaces for GLCM feature computation for microscopic images of microorganisms do not result in an improvement of classification accuracy. Indeed, the best results presented in Table are reached for standard RGB color space, whereas other color channels might decrease classification accuracy from 87.4% to 74.6% using Y’IQ color space and 50 times repeated 10-Fold cross-validation technique and ELM with mish activation function and 1000 neurons on Dataset1. Gray level run length matrix (GLRLM) In general GLCM is not the only method of texture features extraction. Commonly GLRLM is used together with GLCM. GLRLM creates features that describe spatial relationships between pixels having values of a given length of co-occurrence sequences and follows the concept of a run-length encoding outlined below. Non-overlapping sliding of window W similarly to GLCM is performed on the image quantized to the q levels. GLRLM is defined as four matrices [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }$$ each of size [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q r$$ , where r is a maximal run-length of co-occurrence pixels. Rows of that matrix are described as gray level values [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i = 0, ,q-1$$ and columns as run-lengths [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r_j=1, ,r$$ . For every obtained W run-length of co-occurrence pixel values q for direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and maximal run-length r is computed. Assume [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{11}, ,w_{mn}$$ form the coefficients of particular W . The directions are similarly defined as in GLCM method (see Fig. ). In GLRLM the whole window W for each of the gray levels [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i$$ is analyzed. To illustrate the latter for a given W with [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = 45^{ }$$ a maximal run-length [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r = 3$$ and take [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i = 0$$ were set. Then number of run-length sequences of length [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r_j=1, r$$ is noted. Based on such choice, the values in [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{45^{ }}$$ are incremented e.g. there are 2 sequences of length 2 (means that some row in 2nd column should be enlarged) for gray level equal 0 (first row), so coefficient of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{45^{ }}[1,2]$$ is increased by 2 (see Fig. ). Similarly, all other entries of matrix [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^$$ are incremented for every window W based on input image and for every gray level [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i$$ . Consequently all coefficients of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^$$ are divided by qr . As a result matrix [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^$$ contains now probabilities of co-occurrence gray levels [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i$$ for run-length [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r_j$$ and maximum run-length r for a given angle [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ . Based on computed probabilities the following statistical measures for GLRLM (see Zhou et al. ) are calculated: Short Run Emphasis (SRE), Long Run Emphasis, Grey Level Non-uniformity, Run Length Non-uniformity, Run Percentage, Low Grey Level Run Emphasis, High Grey Level Run Emphasis, Short Run Low Gray Level Emphasis, Short Run High Gray Level Emphasis, Long Run Low Gray Level Emphasis and Long Run High Gray Level Emphasis . Similarly to GLCM, each of those features is calculated for four different directions [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and the final value of them is defined as e.g. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$SRE = 0.25(SRE^{0^{ }}+SRE^{45^{ }}+SRE^{90^{ }}+SRE^{135^{ }})$$ . Experiments to improve classification accuracy based on GLRLM features are conducted not to convert image to grayscale, but to calculate GLRLM on each of the color channels separately. Unfortunately, the same approach in case of GLRLM did not improve the results. A comparison is performed between the use of GLRLM feature generation based on images [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ . Surprisingly, better results are achieved using image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ , whereas conducting the same experiment on GLCM do not lead to further improvement in the results. Feature generation summary In total for each dataset the same set of features organized into three color spaces: RGB, I1I2I3, and XYZ (as explained in details in section before). The following abbreviations apply: RGB: Red (R), Green (G), Blue (B) I1I2I3: Intensity channels I1, I2, I3 XYZ: X, Y, Z color space Each color channel has the following statistical features: I1, I2, I3: Mean, Std, Mean (H), Std (H), Kurtosis (H), Skewness (H), Entropy (H), Energy (H) X, Y, Z: Mean, Std, Kurtosis, Skewness GLCM (Gray-Level Co-occurrence Matrix) features: Contrast, Correlation, Energy, Entropy, Homogeneity, Autocorrelation, Cluster Prominence, Inverse Difference, Dissimilarity (for R, G, B channels) GLRLM (Gray-Level Run Length Matrix) features: SRE, LRE, GLN, RLN, RP, LGRE, HGRE, SRLGE, SRHGE, LRLGE, LRHGE Notably the mathematical equations behind these metrics (e.g., Mean, Kurtosis, Skewness, GLCM, GLRLM) can be obtained from corresponding well-known machine vision literature , . Extreme learning machine The ELM network is a fully connected feed-forward neural network (see Huang et. al. ) designed to solve classification problems. This network consists of one hidden-layer with McCulloch-Pitts neurons . The precise optimal number of hidden-layer neurons is each time evaluated empirically as there is no up-front theoretical method determining such value to maximize the classification accuracy. ELM is characterized by a simple topology structure. In network learning, the utilization of the Moore-Penrose pseudoinverse operation generally involves the application of singular value decomposition (SVD). The latter permits to shorten a time-consuming iterative process of adjusting values of network’s weights. Assume that after image segmentation, feature calculation and feature selection dataset consists of N images. It can be described as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = \{(x_i,t_i)\}_{i=1}^{N}$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$x_i$$ is i -th vector of d features and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ is i -th label of class to which selected image belongs. Then input layer of ELM has d neurons, defining also input matrix X introduced in . The number of output layer neurons coincides with the number of image classes M determining the classification task. Thus the output layer is numerically represented by a vector [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y = \{y_i\}_{i=1}^{M}$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_i  {R}$$ . Note that here the maximal value in vector y at index [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_{max}$$ is understood as assigning an input image to the [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_{max}$$ class. In order to test the ELM network, the i -th class label [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ needs to be appropriately formatted and compared to y . A possible approach is to format [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$\{t_{ij}\}_{j=0}^m$$ that has all values set to zero except one element at p -th index that is set to one. Then [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ can be written as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i = [0, ,1, ,0]$$ where p -th element indicates that i -th image affiliates to p -th class of images ( [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_{i_p}=1$$ ). Only the greatest value of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_i$$ at index [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_i^{max}$$ is chosen. Correct classification is observed if and only if [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i[y_i^{max}] = 1$$ for i -th picture. This concept of formatting is called One-Hot Encoding or [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$1-of-K$$ scheme . Let now L be a number of units in single hidden layer network that is chosen empirically. Since ELM is a fully connected network, weights between input layer and hidden one are represented by a matrix W from , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{ij}$$ denotes the weight between i -th neuron in input layer and j -th neuron of hidden layer. The corresponding bias values are assumed to be represented by a vector [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$b = \{b_i\}_{i=1}^{N}$$ . In ELM network both W and b are filled with values that are randomly chosen from the interval [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$[-1,1]$$ using a uniform distribution. Then the outputs of hidden units are calculated as expressed in matrix H (see , where f is any continuous function [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$f:  {R}  {R}$$ . 1 [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ X =  x_{11} &  & x_{1N}\\  &  &  \\ x_{d1} &  & x_{dN} ,  W =  w_{11} &  & w_{1L}\\  &  &  \\ w_{d1} &  & w_{dL} ,  H =  f( _{i=1}^{d}x_{i1}w_{i1} + b_1) &  & f( _{i=1}^{d}x_{i1}w_{iL} + b_1)\\  &  &  \\ f( _{i=1}^{d}x_{iN}w_{i1}+ b_N) &  & f( _{i=1}^{d}x_{iN}w_{iL} + b_N) . $$ Weights [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ between hidden and output layers are still unknown. In order to obtain them, the equitation [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ has to be evaluated. Indeed since H with probability equal to 1 is irreversible and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$\| H -Y \| = 0$$ (see Huang et al. ), the system [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ cannot be directly solved. Instead, one finds the optimal parameter [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = arg\ \| H -T \| ^2 = H^ T$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ is the Moore-Penrose generalized inverse of H (see Rao and Mitra ) and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$T = {t_i}_{i=1}^{N}$$ describing image appurtenance to the given class formed using concept of one hot encoding mentioned in the above paragraph. Moore-Penrose pseudoinverse of a real matrix H must satisfy following conditions: 1. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$HH^ H = H,$$ 2. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^ HH^ = H^ ,$$ 3. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$(HH^ H)^T = HH^ ,$$ 4. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$(H^ H)^T = H^ H.$$ Pseudoinverse is known to be unique and pseudoinverse of a non-singular quadratic matrix is identical to the ordinary inverse. For matrices that are singular pseudoinverse solves [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y = H$$ in terms of least-squares problem. Calculating [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ gives solution [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ so that [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H$$ is as close as possible to Y in terms of Euclidean distance. Note here that H generically corresponds to either overdetermined or underdetermined system. Additionally, there can be many methods to calculate [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ , one of them recalled here is orthogonal projection method (see Rao and Mitra ). If [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$L>d$$ (number of hidden units is greater than number of input units) then system [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ is underdetermined and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^ = H^T(HH^T)^{-1}$$ , otherwise system [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ is overdetermined and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^ = (H^TH)^{-1}H^T$$ . In practical applications of ELM [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$L d$$ eventuates. In classification tasks L usually totals a few thousands. Then, operation of inverse [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$HH^T$$ could be time-consuming especially taking into account the fact that operation of matrix pseudoinverse cannot be parallelized. To improve learning speed of network QR matrix decomposition is used (see Strang ). QR decomposition of matrix A is defined as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$A=QR$$ , where R is upper triangular and Q is orthogonal ( [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$QQ^T =  {I}$$ ). Using QR decomposition, the matrix [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ can be computed according to: [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ H^{ } = H^T(HH^T)^{-1} = R^TQ^T(R^TQ^TQR)^{-1} = R^TQ^T(R^TR)^{-1} = Q^TR^{-1}. $$ The ELM network gives a solution that minimizes mean square error between predicted [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and actual Y . Thanks to assigning random values to weights and bias between input and hidden layer, the trained network is very resistant to overtraining. The computed solution is not a local minimum of optimization criterion like in MLP with backpropagation algorithm based on gradient descent , . The network training process is much faster than other Machine Learning methods like MLP, RF and SVM. The main drawback of this process relies on finding L and f that maximize classification accuracy.   Microorganisms are grown in Potato Dextrose agar, specifically Merck 1.10130.0500 at 26 degrees Celsius. In the case of Fusarium , Verticillium and Phytophthora the microorganisms are incubated for 7–10 days, while for Trichoderma strain it takes about 4-5 days until conidial spores are produced. In total five different datasets were formed that can be obtained from the corresponding author upon the reasonable request. Dataset1 The initial experiments utilized a dataset comprising 128 images of microorganisms from the genera Verticillium , Trichoderma , Phytophthora and Fusarium , referred to as Dataset1 throughout this manuscript with exemplary images presented in Fig. . This foundational dataset is based on the work of Kruk et al. , which demonstrated the potential for classifying soil-dwelling microorganisms using ML techniques. The primary objective with Dataset1 was to develop a robust image preprocessing and segmentation pipeline that could accurately differentiate the microorganism regions from the background as this is recognized as a crucial step for ultimate model’s performance. Despite being well-suited for preliminary investigations into image segmentation and feature extraction, Dataset1 was relatively small, which limited its ability to ensure the robustness of the solution across diverse scenarios. To overcome these limitations and to provide a comprehensive evaluation of the methodology’s performance, the additional, substantially larger datasets are provided. These new datasets were crucial in validating the generalizability and effectiveness of the proposed image preprocessing, masking and classification methods at scale. Dataset1 was carefully curated by professional microbiologists who were responsible for both sample preparation and image acquisition through microscopy. The primary goal was to generate high-quality images that clearly showcased the morphological characteristics essential for accurate microorganism identification. To achieve this, each image was meticulously focused to ensure that the unique features of each genus were prominently visible. The process began with the cultivation of the microorganisms, followed by the careful transfer biological material onto microscope slides. To facilitate optimal imaging, a drop of distilled water was added to dilute the biological material and a cover slip was placed on top (see Fig. ). These prepared slides were then subjected to microscopic examination, allowing the microbiologists to capture detailed images at different magnification levels that would serve as the basis for the ML analysis. Given that the model in question solely relies on image data for input without the incorporation of additional multimodal features, the quality and clarity of these images were paramount. However, this approach introduces a degree of bias, as image acquisition process rely heavily on the expertise of the microbiologists. This necessitates the involvement of experts throughout the pipeline to ensure that the images accurately represent the microorganisms’ traits for effective model training and evaluation. Dataset2 Dataset2 builds upon the foundation of Dataset1 by expanding its scope and content. This expanded dataset includes additional images of the original genera ( Trichoderma and Fusarium ) and introduces images of a new genus, Purpureocillium (see Fig. ). Prepared using the same procedures as Dataset1, Dataset2 effectively consolidates and enhances the original set, resulting in a comprehensive collection of 304 images. Dataset3 and Dataset4 Upon evaluating Dataset1 the accuracy of the generated segmentation masks is assessed for the suitability of the calculated features and the overall performance of the applied ML methods. This comprehensive analysis highlighted the need for a larger, more automated dataset to minimize bias and enhance scalability. Thus, the third dataset, referred to in this paper as Dataset3 was prepared, aiming to automate the entire image acquisition process and build a significantly larger dataset without the direct involvement of experts. The critical point in creating Dataset3 was the elimination of expert bias inherent in manually capturing images, which, while valuable for understanding morphological traits, does not reflect real-world scenarios where a fully automated system is the target. The approach involved using an automated microscope equipped with a movable mechanical stage controlled by built-in servo motors and software drivers. This setup allowed the microscope to systematically capture images across a predefined area without overlap. During initial trials, the challenges with out-of-focus images due to the varying depths (Z-axis) of microorganisms at high magnification levels and slight imperfections on the slide surfaces were encountered. Although, the autofocus feature helped address this issue, it occasionally failed to achieve perfect focus, leading to some variability in image quality. Despite these challenges, the automated system substantially reduced the reliance on expert involvement and increased the efficiency of image acquisition. To ensure comparability, the samples for Dataset3 were cultivated similarly to those in Dataset1. Two methods of material transfer onto microscope slides were examined. The first method, consistent with Dataset1, involved placing the biological material on a slide, adding a drop of water and covering it with a slip. This approach created Dataset3. The second method involved smearing the material directly onto the slide without a cover slip, which formed a new dataset, referred to as Dataset4 (see Fig. b). Both Dataset3 and Dataset4 were collected without requiring biological expertise for each image. Instead, a trained man handled the sample preparation and setup of the microscope, which included positioning the sample on the mechanical stage, using software to target areas of interest and initiating the automated imaging process with autofocus. This approach simulates a real-world scenario where minimal expert intervention is needed. Sample images from both datasets are presented in Fig. . Looking forward, there is potential for further automation, such as using robots for material transfer and advanced software for precise targeting, which would enhance the efficiency and reliability of the image acquisition process. These advancements are crucial for the ongoing development and scalability of the system. Dataset5 The fifth dataset, referred to as Dataset5, was generated following a procedure similar to that of Dataset4 (Fig. ). For Dataset4, images were captured in a non-overlapping, systematic manner to cover the entire sample area, as illustrated in Fig. b. Each image was taken using the microscope’s autofocus feature, which aimed to optimize focus for each shot. However, while autofocus provided a reasonable solution, it occasionally failed to achieve perfect focus across all images. To address the limitations of autofocus, Dataset5 incorporated a depth imaging technique. This method involved capturing a series of images at varying focal distances (the distance between the sample and the camera lens) for each section of the sample. These multiple images, each focused at a different depth, were then automatically stacked to produce a single composite image. This stacking process significantly enhances the overall image quality by ensuring that all areas of the sample are in focus, thus mitigating the problem of focal inconsistencies observed with single-focus images. Although this depth imaging approach enhances image clarity significantly, it requires a slower data acquisition process due to the need to capture and stack multiple images at different focal depths for each area of the sample. For Dataset5, the additional challenges were encountered when attempting to use the depth imaging technique in conjunction with the traditional method of adding a drop of water between the microscope slide and the coverslip that was applied for Dataset3. This method, commonly used to stabilize the sample for conventional imaging, proved unsuitable for depth imaging. The water layer caused instability, leading to slight movements of fungal elements within the sample. Consequently, these movements resulted in blurred or misaligned images when stacking the focal series, thereby undermining the intended clarity of the depth imaging process. Overview of the datasets Initially, Dataset1 was prepared for preliminary research, focusing on data segmentation and feature calculation. This dataset included 128 images from four genera: Verticillium , Trichoderma , Phytophthora and Fusarium . To test and to enhance the robustness of the solution, this dataset was expanded by creating Dataset2. This expanded dataset added more images for Trichoderma and Fusarium and incorporated a new genus, Purpureocillium , increasing the total to 303 images. Subsequently, with the implementation of an automated image acquisition system, three additional datasets were developed: Dataset3, Dataset4 and Dataset5. These datasets were captured using an automated microscope, allowing to gather a significantly larger and more diverse collection of images. Unfortunately, Phytophthora could not be included in the automated datasets due to deterioration of all samples held prior to the establishment of the automated imaging system. Collectively, these efforts resulted in the creation of the largest database of soil microorganism images captured through microscopy for the ML purposes, setting a new standard in the literature. Datasets with all details are provided in Table . The findings presented in this research are crucial due to the successful replication of the entire machine vision algorithm and ML methods across a range of datasets. These datasets vary significantly in size and were acquired using both automated and manual methods, as well as through two distinct types of microscopes. Such replication demonstrates that the system retains its performance consistently across different conditions, indicating a strong resistance to overfitting. For statistical proofs of this differences please refer to section Feature analysis . This robustness is evident from the ability to maintain high accuracy and reliability in image analysis and classification, regardless of how the data was collected. The diverse operational settings in which the datasets were obtained, including automated and non-automated means, further validate the adaptability of the system. This versatility ensures that the system may be successfully deployed in a variety of real-world circumstances where consistency and reliability are critical. By proving effective across a broad spectrum of datasets and acquisition techniques, the system showcases its potential for scalability and generalization in the field of microbial image analysis.   The initial experiments utilized a dataset comprising 128 images of microorganisms from the genera Verticillium , Trichoderma , Phytophthora and Fusarium , referred to as Dataset1 throughout this manuscript with exemplary images presented in Fig. . This foundational dataset is based on the work of Kruk et al. , which demonstrated the potential for classifying soil-dwelling microorganisms using ML techniques. The primary objective with Dataset1 was to develop a robust image preprocessing and segmentation pipeline that could accurately differentiate the microorganism regions from the background as this is recognized as a crucial step for ultimate model’s performance. Despite being well-suited for preliminary investigations into image segmentation and feature extraction, Dataset1 was relatively small, which limited its ability to ensure the robustness of the solution across diverse scenarios. To overcome these limitations and to provide a comprehensive evaluation of the methodology’s performance, the additional, substantially larger datasets are provided. These new datasets were crucial in validating the generalizability and effectiveness of the proposed image preprocessing, masking and classification methods at scale. Dataset1 was carefully curated by professional microbiologists who were responsible for both sample preparation and image acquisition through microscopy. The primary goal was to generate high-quality images that clearly showcased the morphological characteristics essential for accurate microorganism identification. To achieve this, each image was meticulously focused to ensure that the unique features of each genus were prominently visible. The process began with the cultivation of the microorganisms, followed by the careful transfer biological material onto microscope slides. To facilitate optimal imaging, a drop of distilled water was added to dilute the biological material and a cover slip was placed on top (see Fig. ). These prepared slides were then subjected to microscopic examination, allowing the microbiologists to capture detailed images at different magnification levels that would serve as the basis for the ML analysis. Given that the model in question solely relies on image data for input without the incorporation of additional multimodal features, the quality and clarity of these images were paramount. However, this approach introduces a degree of bias, as image acquisition process rely heavily on the expertise of the microbiologists. This necessitates the involvement of experts throughout the pipeline to ensure that the images accurately represent the microorganisms’ traits for effective model training and evaluation.   Dataset2 builds upon the foundation of Dataset1 by expanding its scope and content. This expanded dataset includes additional images of the original genera ( Trichoderma and Fusarium ) and introduces images of a new genus, Purpureocillium (see Fig. ). Prepared using the same procedures as Dataset1, Dataset2 effectively consolidates and enhances the original set, resulting in a comprehensive collection of 304 images.   Upon evaluating Dataset1 the accuracy of the generated segmentation masks is assessed for the suitability of the calculated features and the overall performance of the applied ML methods. This comprehensive analysis highlighted the need for a larger, more automated dataset to minimize bias and enhance scalability. Thus, the third dataset, referred to in this paper as Dataset3 was prepared, aiming to automate the entire image acquisition process and build a significantly larger dataset without the direct involvement of experts. The critical point in creating Dataset3 was the elimination of expert bias inherent in manually capturing images, which, while valuable for understanding morphological traits, does not reflect real-world scenarios where a fully automated system is the target. The approach involved using an automated microscope equipped with a movable mechanical stage controlled by built-in servo motors and software drivers. This setup allowed the microscope to systematically capture images across a predefined area without overlap. During initial trials, the challenges with out-of-focus images due to the varying depths (Z-axis) of microorganisms at high magnification levels and slight imperfections on the slide surfaces were encountered. Although, the autofocus feature helped address this issue, it occasionally failed to achieve perfect focus, leading to some variability in image quality. Despite these challenges, the automated system substantially reduced the reliance on expert involvement and increased the efficiency of image acquisition. To ensure comparability, the samples for Dataset3 were cultivated similarly to those in Dataset1. Two methods of material transfer onto microscope slides were examined. The first method, consistent with Dataset1, involved placing the biological material on a slide, adding a drop of water and covering it with a slip. This approach created Dataset3. The second method involved smearing the material directly onto the slide without a cover slip, which formed a new dataset, referred to as Dataset4 (see Fig. b). Both Dataset3 and Dataset4 were collected without requiring biological expertise for each image. Instead, a trained man handled the sample preparation and setup of the microscope, which included positioning the sample on the mechanical stage, using software to target areas of interest and initiating the automated imaging process with autofocus. This approach simulates a real-world scenario where minimal expert intervention is needed. Sample images from both datasets are presented in Fig. . Looking forward, there is potential for further automation, such as using robots for material transfer and advanced software for precise targeting, which would enhance the efficiency and reliability of the image acquisition process. These advancements are crucial for the ongoing development and scalability of the system.   The fifth dataset, referred to as Dataset5, was generated following a procedure similar to that of Dataset4 (Fig. ). For Dataset4, images were captured in a non-overlapping, systematic manner to cover the entire sample area, as illustrated in Fig. b. Each image was taken using the microscope’s autofocus feature, which aimed to optimize focus for each shot. However, while autofocus provided a reasonable solution, it occasionally failed to achieve perfect focus across all images. To address the limitations of autofocus, Dataset5 incorporated a depth imaging technique. This method involved capturing a series of images at varying focal distances (the distance between the sample and the camera lens) for each section of the sample. These multiple images, each focused at a different depth, were then automatically stacked to produce a single composite image. This stacking process significantly enhances the overall image quality by ensuring that all areas of the sample are in focus, thus mitigating the problem of focal inconsistencies observed with single-focus images. Although this depth imaging approach enhances image clarity significantly, it requires a slower data acquisition process due to the need to capture and stack multiple images at different focal depths for each area of the sample. For Dataset5, the additional challenges were encountered when attempting to use the depth imaging technique in conjunction with the traditional method of adding a drop of water between the microscope slide and the coverslip that was applied for Dataset3. This method, commonly used to stabilize the sample for conventional imaging, proved unsuitable for depth imaging. The water layer caused instability, leading to slight movements of fungal elements within the sample. Consequently, these movements resulted in blurred or misaligned images when stacking the focal series, thereby undermining the intended clarity of the depth imaging process.   Initially, Dataset1 was prepared for preliminary research, focusing on data segmentation and feature calculation. This dataset included 128 images from four genera: Verticillium , Trichoderma , Phytophthora and Fusarium . To test and to enhance the robustness of the solution, this dataset was expanded by creating Dataset2. This expanded dataset added more images for Trichoderma and Fusarium and incorporated a new genus, Purpureocillium , increasing the total to 303 images. Subsequently, with the implementation of an automated image acquisition system, three additional datasets were developed: Dataset3, Dataset4 and Dataset5. These datasets were captured using an automated microscope, allowing to gather a significantly larger and more diverse collection of images. Unfortunately, Phytophthora could not be included in the automated datasets due to deterioration of all samples held prior to the establishment of the automated imaging system. Collectively, these efforts resulted in the creation of the largest database of soil microorganism images captured through microscopy for the ML purposes, setting a new standard in the literature. Datasets with all details are provided in Table . The findings presented in this research are crucial due to the successful replication of the entire machine vision algorithm and ML methods across a range of datasets. These datasets vary significantly in size and were acquired using both automated and manual methods, as well as through two distinct types of microscopes. Such replication demonstrates that the system retains its performance consistently across different conditions, indicating a strong resistance to overfitting. For statistical proofs of this differences please refer to section Feature analysis . This robustness is evident from the ability to maintain high accuracy and reliability in image analysis and classification, regardless of how the data was collected. The diverse operational settings in which the datasets were obtained, including automated and non-automated means, further validate the adaptability of the system. This versatility ensures that the system may be successfully deployed in a variety of real-world circumstances where consistency and reliability are critical. By proving effective across a broad spectrum of datasets and acquisition techniques, the system showcases its potential for scalability and generalization in the field of microbial image analysis.   All input images are automatically preprocessed in order to improve the segmentation phase which in turn has an essential impact on final classification accuracy. Indeed, the primary challenge associated with microscopic images lies in managing noise interference and low contrast impact. Another significant consideration is the presence of scale labels positioned at the bottom-left corner of each image for Dataset1 and Dataset2 only. It is important to note that the colors of these labels may vary which potentially brings distortion to the generated results. The latter occurs due to the automatic color selection of the labels by a microscopic software. Image masks that differentiate image to label image foreground and background are prepared by thresholding, whose values are empirically chosen. Consequently, the labels are concealed by replacing their area with the average color value of the complete image. The following stage involves converting the given image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}$$ in RGB color space to grayscale forming image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{gray}$$ , which can be subjected to further processing. To mitigate image noise, this study applies Gaussian and Wiener filters . Finally, to improve the clarity of image features and to aid in distinguishing microorganisms from the background, the Contrast Limited Adaptive Histogram Equalization (CLAHE) algorithm is employed.   Majority of the microscopic images in the datasets exhibit an imbalanced ratio between the area of the background and the area of the region of interest (ROI) that contains the soil microorganisms. Most of the area on microscopic pictures is covered by background information irrelevant to the classification task e.g. test-tube, medium for the cultivation of bacteria or minor impurities. Therefore, background has to be separated from the ROI. To address this issue, a binary image mask is created, assigning a value of 1 to each pixel that belongs to the ROI and a value of 0 to each pixel that belongs to the background. One of the most commonly used methods to perform image segmentation is an Otsu method that calculates one threshold value globally for the whole image upon minimizing weighted sum of variances within ROI and background pixel values and maximizing weighted sum of variances between them. According to , the Otsu method is known to produce weak results when the area of the ROI is imbalanced in comparison to the background area. Due to the above limitation a method called Adaptive Image Thresholding is used. The mask obtained using adaptive image threshold needs further image processing using morphological operations. Open operation is used to remove small objects in image and close operation is applied to join together objects that are close to each other. Both operations are applied to remove segmentation imperfections caused by the noise occurrence within original image (the so-called salt and pepper noise removal—see ). Despite incorporating morphological operations, the thresholding results remain unsatisfactory. Classification accuracy of the ELM classifier applied on features computed for the generated mask renders 82.7% for Dataset1. Unsatisfactory results are caused by the fact that input images are very complex as having different structures and textures. The segmentation accuracy is enhanced up to 87.4%, thus improving the performance of the classifier upon applying k -means algorithm . Mask M is multiplied using Hadamard product with [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ forming combined image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ . On that image k -means method with [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$k=3$$ is run. As a result the original image is separated into three classes: i.e. a background and ROI area subsequently divided to two disjoint subclasses. Pixels that are classified to the brightest centroid represent a new, smaller area of ROI. The pixels with the highest contrast are assumed to have the most relevant information to distinguish microorganisms. In contrast, ROI pixels assigned to the centroid with lower values are classified in fact to the background as it may contain data disrupting the final results. In this manner a new mask M is generated. Finally, on a given mask, open and close morphological operations are applied. All crucial image segmentation stjpg contributing to the above algorithm are illustrated in Fig. and the results from various datasets in Fig. .   Image features are computed as descriptors to measure statistical attributes of an image. The primary objective is to devise features that can effectively discriminate the different genera of microorganisms, enabling the formation of a robust classifier. Thus, color and texture-based handcrafted features are calculated in this study that are explained in the forthcoming subsections. Color features Color features are computed on combined images [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ generated by Hadamard product of original image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ with its mask M introduced in the previous section. Conversion of an image from RGB color values to a different color space prior to feature calculation can enhance classification accuracy . Indeed, experiments show such increase once conversion from RGB to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ color space is applied. First set of color features contains 8 statistical features, each calculated separately for three image color channels of color space [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ to which [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ was previously converted. Every image pixel of coordinates x and y in RGB image is a vector containing 3 values for the respective RGB color space i.e. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_{RGB}[x,y] = [R[x,y],G[x,y],B[x,y]]$$ . Conversion from RGB to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ is made upon applying the following operation: [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ I_{I_1I_2I_3}[x,y] = I_{RGB}[x,y] 1/3 & 1/3 & 1/3\\ 1/2 & 0 & -1/2\\ -1/4 & 1/2 & -1/4 . $$ Histograms for images [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ (denoted by H and h , respectively, generated upon conversion to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ color space) are calculated for every color space channel. Image histogram is a vector h , where h [ i ] is a number of occurrences of image pixels with the corresponding intensity value ranging within [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0 i 255$$ . Having both histograms the following features are calculated: mean and standard deviation of H; mean, standard deviation, kurtosis, skewness, entropy and energy of h . The second set of color features is based on image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ converted into XYZ color space. In contrast to the previous case, the calculated measures rely directly on image pixels, not histogram distribution. The latter permits to extract additional information on pixel spatial relationships which in turn is unavailable from image histogram . The XYZ color space is selected after comparison of classification accuracy for the different color spaces. For every X , Y and Z color channel four statistical measures are calculated separately i.e. mean, standard deviation, kurtosis and skewness yielding in total 12 features.   Color features are computed on combined images [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ generated by Hadamard product of original image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ with its mask M introduced in the previous section. Conversion of an image from RGB color values to a different color space prior to feature calculation can enhance classification accuracy . Indeed, experiments show such increase once conversion from RGB to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ color space is applied. First set of color features contains 8 statistical features, each calculated separately for three image color channels of color space [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ to which [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ was previously converted. Every image pixel of coordinates x and y in RGB image is a vector containing 3 values for the respective RGB color space i.e. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_{RGB}[x,y] = [R[x,y],G[x,y],B[x,y]]$$ . Conversion from RGB to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ is made upon applying the following operation: [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ I_{I_1I_2I_3}[x,y] = I_{RGB}[x,y] 1/3 & 1/3 & 1/3\\ 1/2 & 0 & -1/2\\ -1/4 & 1/2 & -1/4 . $$ Histograms for images [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ (denoted by H and h , respectively, generated upon conversion to [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_1I_2I_3$$ color space) are calculated for every color space channel. Image histogram is a vector h , where h [ i ] is a number of occurrences of image pixels with the corresponding intensity value ranging within [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0 i 255$$ . Having both histograms the following features are calculated: mean and standard deviation of H; mean, standard deviation, kurtosis, skewness, entropy and energy of h . The second set of color features is based on image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ converted into XYZ color space. In contrast to the previous case, the calculated measures rely directly on image pixels, not histogram distribution. The latter permits to extract additional information on pixel spatial relationships which in turn is unavailable from image histogram . The XYZ color space is selected after comparison of classification accuracy for the different color spaces. For every X , Y and Z color channel four statistical measures are calculated separately i.e. mean, standard deviation, kurtosis and skewness yielding in total 12 features.   Image texture defines spatial relations of chrominance or luminance of the image pixels. Grey-level co-occurrence matrix (GLCM) GLCM introduced by Haralick is a method to extract spatial relations of pixels for a given angle direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and maximal distance d . To apply the GLCM algorithm, the first step involves quantizing the image into q levels, resulting in the creation of an image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q$$ . Next, four co-occurrence matrices [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }_{q q}$$ are defined, each computed for a different angle [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ (i.e., 0, 45, 90 degrees and 135 degrees). Here all coefficients of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }_{q q}$$ are initialized to zero. Examination of pixels relations is performed over their local neighborhood using window W of size [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$m n$$ pixels. In sequel, W is iteratively moved throughout the image I assuming that coordinates of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q$$ are positioned in the top left corner. (see Fig. ). At the first iteration of algorithm matrix W contains pixels of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q[y,x]$$ with indices satisfying [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<y m$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<x n$$ . Next W is moved horizontally so that it does not overlap with the previous window - e.g. indices of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q[y,x]$$ for the second window range over [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<y m$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$n<x  2n$$ . If such left-to-right shift of W is not anymore possible such window is shifted vertically to underneath-left boundary position (e.g. after the first horizontal shifting movement the indices run over [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$m<y 2m$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<x n$$ as illustrated in Fig. (see window III). GLCM algorithm terminates once none of the above stjpg is further possible (i.e. W is too close to the right bottom corner of the I ). Every obtained W is the subject of further calculations. Co-occurrences of pixel values [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i=0, , q-1$$ are noted for a direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = 0^{ }, 45^{ }, 90^{ }, 135^{ }$$ and maximal distance d . Let [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{11}, ,w_{mn}$$ be coefficients of a given matrix W . Iterating through all of the coefficients one analyzes spatial correlations in their neighborhood. Currently analyzed pixel (coefficient of W in certain iteration) is named pixel of interest p . Definition of directions for maximal distance [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$d=2$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$p = w_{42}$$ is presented in Fig. . Assume p is [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{42}=0$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ =0^{ }$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$d=2$$ . Rows in [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ store the value of p and columns values of co-occurrence pixels being in relation to p . Since [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{42}=0$$ , the only coefficients in the first row of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ can be incremented. Then, based on Fig. , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{43}$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{44}$$ should be now analyzed. Coefficient [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{43}=0$$ implies that value in [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ in the first column ( [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}[1,1]$$ ) is incremented. Analogously, for [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{43}=3$$ [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}[1,4]$$ is increased. Both stjpg are presented in Fig. . In the same manner calculations for every W are performed analyzing all of the [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{ij}$$ for a direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ . Algorithm’s result is four matrices: [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{45^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{90^{ }}$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{135^{ }}$$ whose coefficients are then converted to the probabilities by dividing each of its values by mn . In the next step, the following statistical measures called GLCM features based on coefficients in matrices [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }$$ can be computed: Contrast, Correlation, Energy, Homogeneity, Autocorrelation, Cluster Prominence, Inverse Difference and Dissimilarity . Noticeably each of these features is calculable separately for each [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = 0^{ }, 45^{ }, 90^{ }, 135^{ }$$ . Having computed [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{0^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{45^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{90^{ }}$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{135^{ }}$$ their mean [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy = 0.25(Energy^{0^{ }}+Energy^{45^{ }}+Energy^{90^{ }}+Energy^{135^{ }})$$ is next derived. Originally GLCM based features are calculated on the grayscale image i.e. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_t^{gray}$$ . A further improvement is obtained here once GLCM features are computed separately for each of the RGB image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_t^{RGB}$$ channel yielding in total 24 features. In contrast to the work by Duan et al. , which proposes an extension that utilizes images in a different color space than RGB for GLCM feature computation, the experiments incorporating different color spaces for GLCM feature computation for microscopic images of microorganisms do not result in an improvement of classification accuracy. Indeed, the best results presented in Table are reached for standard RGB color space, whereas other color channels might decrease classification accuracy from 87.4% to 74.6% using Y’IQ color space and 50 times repeated 10-Fold cross-validation technique and ELM with mish activation function and 1000 neurons on Dataset1. Gray level run length matrix (GLRLM) In general GLCM is not the only method of texture features extraction. Commonly GLRLM is used together with GLCM. GLRLM creates features that describe spatial relationships between pixels having values of a given length of co-occurrence sequences and follows the concept of a run-length encoding outlined below. Non-overlapping sliding of window W similarly to GLCM is performed on the image quantized to the q levels. GLRLM is defined as four matrices [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }$$ each of size [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q r$$ , where r is a maximal run-length of co-occurrence pixels. Rows of that matrix are described as gray level values [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i = 0, ,q-1$$ and columns as run-lengths [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r_j=1, ,r$$ . For every obtained W run-length of co-occurrence pixel values q for direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and maximal run-length r is computed. Assume [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{11}, ,w_{mn}$$ form the coefficients of particular W . The directions are similarly defined as in GLCM method (see Fig. ). In GLRLM the whole window W for each of the gray levels [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i$$ is analyzed. To illustrate the latter for a given W with [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = 45^{ }$$ a maximal run-length [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r = 3$$ and take [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i = 0$$ were set. Then number of run-length sequences of length [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r_j=1, r$$ is noted. Based on such choice, the values in [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{45^{ }}$$ are incremented e.g. there are 2 sequences of length 2 (means that some row in 2nd column should be enlarged) for gray level equal 0 (first row), so coefficient of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{45^{ }}[1,2]$$ is increased by 2 (see Fig. ). Similarly, all other entries of matrix [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^$$ are incremented for every window W based on input image and for every gray level [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i$$ . Consequently all coefficients of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^$$ are divided by qr . As a result matrix [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^$$ contains now probabilities of co-occurrence gray levels [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i$$ for run-length [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r_j$$ and maximum run-length r for a given angle [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ . Based on computed probabilities the following statistical measures for GLRLM (see Zhou et al. ) are calculated: Short Run Emphasis (SRE), Long Run Emphasis, Grey Level Non-uniformity, Run Length Non-uniformity, Run Percentage, Low Grey Level Run Emphasis, High Grey Level Run Emphasis, Short Run Low Gray Level Emphasis, Short Run High Gray Level Emphasis, Long Run Low Gray Level Emphasis and Long Run High Gray Level Emphasis . Similarly to GLCM, each of those features is calculated for four different directions [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and the final value of them is defined as e.g. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$SRE = 0.25(SRE^{0^{ }}+SRE^{45^{ }}+SRE^{90^{ }}+SRE^{135^{ }})$$ . Experiments to improve classification accuracy based on GLRLM features are conducted not to convert image to grayscale, but to calculate GLRLM on each of the color channels separately. Unfortunately, the same approach in case of GLRLM did not improve the results. A comparison is performed between the use of GLRLM feature generation based on images [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ . Surprisingly, better results are achieved using image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ , whereas conducting the same experiment on GLCM do not lead to further improvement in the results. Feature generation summary In total for each dataset the same set of features organized into three color spaces: RGB, I1I2I3, and XYZ (as explained in details in section before). The following abbreviations apply: RGB: Red (R), Green (G), Blue (B) I1I2I3: Intensity channels I1, I2, I3 XYZ: X, Y, Z color space Each color channel has the following statistical features: I1, I2, I3: Mean, Std, Mean (H), Std (H), Kurtosis (H), Skewness (H), Entropy (H), Energy (H) X, Y, Z: Mean, Std, Kurtosis, Skewness GLCM (Gray-Level Co-occurrence Matrix) features: Contrast, Correlation, Energy, Entropy, Homogeneity, Autocorrelation, Cluster Prominence, Inverse Difference, Dissimilarity (for R, G, B channels) GLRLM (Gray-Level Run Length Matrix) features: SRE, LRE, GLN, RLN, RP, LGRE, HGRE, SRLGE, SRHGE, LRLGE, LRHGE Notably the mathematical equations behind these metrics (e.g., Mean, Kurtosis, Skewness, GLCM, GLRLM) can be obtained from corresponding well-known machine vision literature , . Extreme learning machine The ELM network is a fully connected feed-forward neural network (see Huang et. al. ) designed to solve classification problems. This network consists of one hidden-layer with McCulloch-Pitts neurons . The precise optimal number of hidden-layer neurons is each time evaluated empirically as there is no up-front theoretical method determining such value to maximize the classification accuracy. ELM is characterized by a simple topology structure. In network learning, the utilization of the Moore-Penrose pseudoinverse operation generally involves the application of singular value decomposition (SVD). The latter permits to shorten a time-consuming iterative process of adjusting values of network’s weights. Assume that after image segmentation, feature calculation and feature selection dataset consists of N images. It can be described as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = \{(x_i,t_i)\}_{i=1}^{N}$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$x_i$$ is i -th vector of d features and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ is i -th label of class to which selected image belongs. Then input layer of ELM has d neurons, defining also input matrix X introduced in . The number of output layer neurons coincides with the number of image classes M determining the classification task. Thus the output layer is numerically represented by a vector [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y = \{y_i\}_{i=1}^{M}$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_i  {R}$$ . Note that here the maximal value in vector y at index [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_{max}$$ is understood as assigning an input image to the [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_{max}$$ class. In order to test the ELM network, the i -th class label [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ needs to be appropriately formatted and compared to y . A possible approach is to format [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$\{t_{ij}\}_{j=0}^m$$ that has all values set to zero except one element at p -th index that is set to one. Then [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ can be written as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i = [0, ,1, ,0]$$ where p -th element indicates that i -th image affiliates to p -th class of images ( [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_{i_p}=1$$ ). Only the greatest value of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_i$$ at index [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_i^{max}$$ is chosen. Correct classification is observed if and only if [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i[y_i^{max}] = 1$$ for i -th picture. This concept of formatting is called One-Hot Encoding or [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$1-of-K$$ scheme . Let now L be a number of units in single hidden layer network that is chosen empirically. Since ELM is a fully connected network, weights between input layer and hidden one are represented by a matrix W from , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{ij}$$ denotes the weight between i -th neuron in input layer and j -th neuron of hidden layer. The corresponding bias values are assumed to be represented by a vector [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$b = \{b_i\}_{i=1}^{N}$$ . In ELM network both W and b are filled with values that are randomly chosen from the interval [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$[-1,1]$$ using a uniform distribution. Then the outputs of hidden units are calculated as expressed in matrix H (see , where f is any continuous function [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$f:  {R}  {R}$$ . 1 [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ X =  x_{11} &  & x_{1N}\\  &  &  \\ x_{d1} &  & x_{dN} ,  W =  w_{11} &  & w_{1L}\\  &  &  \\ w_{d1} &  & w_{dL} ,  H =  f( _{i=1}^{d}x_{i1}w_{i1} + b_1) &  & f( _{i=1}^{d}x_{i1}w_{iL} + b_1)\\  &  &  \\ f( _{i=1}^{d}x_{iN}w_{i1}+ b_N) &  & f( _{i=1}^{d}x_{iN}w_{iL} + b_N) . $$ Weights [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ between hidden and output layers are still unknown. In order to obtain them, the equitation [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ has to be evaluated. Indeed since H with probability equal to 1 is irreversible and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$\| H -Y \| = 0$$ (see Huang et al. ), the system [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ cannot be directly solved. Instead, one finds the optimal parameter [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = arg\ \| H -T \| ^2 = H^ T$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ is the Moore-Penrose generalized inverse of H (see Rao and Mitra ) and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$T = {t_i}_{i=1}^{N}$$ describing image appurtenance to the given class formed using concept of one hot encoding mentioned in the above paragraph. Moore-Penrose pseudoinverse of a real matrix H must satisfy following conditions: 1. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$HH^ H = H,$$ 2. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^ HH^ = H^ ,$$ 3. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$(HH^ H)^T = HH^ ,$$ 4. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$(H^ H)^T = H^ H.$$ Pseudoinverse is known to be unique and pseudoinverse of a non-singular quadratic matrix is identical to the ordinary inverse. For matrices that are singular pseudoinverse solves [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y = H$$ in terms of least-squares problem. Calculating [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ gives solution [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ so that [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H$$ is as close as possible to Y in terms of Euclidean distance. Note here that H generically corresponds to either overdetermined or underdetermined system. Additionally, there can be many methods to calculate [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ , one of them recalled here is orthogonal projection method (see Rao and Mitra ). If [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$L>d$$ (number of hidden units is greater than number of input units) then system [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ is underdetermined and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^ = H^T(HH^T)^{-1}$$ , otherwise system [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ is overdetermined and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^ = (H^TH)^{-1}H^T$$ . In practical applications of ELM [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$L d$$ eventuates. In classification tasks L usually totals a few thousands. Then, operation of inverse [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$HH^T$$ could be time-consuming especially taking into account the fact that operation of matrix pseudoinverse cannot be parallelized. To improve learning speed of network QR matrix decomposition is used (see Strang ). QR decomposition of matrix A is defined as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$A=QR$$ , where R is upper triangular and Q is orthogonal ( [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$QQ^T =  {I}$$ ). Using QR decomposition, the matrix [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ can be computed according to: [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ H^{ } = H^T(HH^T)^{-1} = R^TQ^T(R^TQ^TQR)^{-1} = R^TQ^T(R^TR)^{-1} = Q^TR^{-1}. $$ The ELM network gives a solution that minimizes mean square error between predicted [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and actual Y . Thanks to assigning random values to weights and bias between input and hidden layer, the trained network is very resistant to overtraining. The computed solution is not a local minimum of optimization criterion like in MLP with backpropagation algorithm based on gradient descent , . The network training process is much faster than other Machine Learning methods like MLP, RF and SVM. The main drawback of this process relies on finding L and f that maximize classification accuracy.   GLCM introduced by Haralick is a method to extract spatial relations of pixels for a given angle direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and maximal distance d . To apply the GLCM algorithm, the first step involves quantizing the image into q levels, resulting in the creation of an image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q$$ . Next, four co-occurrence matrices [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }_{q q}$$ are defined, each computed for a different angle [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ (i.e., 0, 45, 90 degrees and 135 degrees). Here all coefficients of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }_{q q}$$ are initialized to zero. Examination of pixels relations is performed over their local neighborhood using window W of size [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$m n$$ pixels. In sequel, W is iteratively moved throughout the image I assuming that coordinates of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q$$ are positioned in the top left corner. (see Fig. ). At the first iteration of algorithm matrix W contains pixels of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q[y,x]$$ with indices satisfying [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<y m$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<x n$$ . Next W is moved horizontally so that it does not overlap with the previous window - e.g. indices of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_q[y,x]$$ for the second window range over [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<y m$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$n<x  2n$$ . If such left-to-right shift of W is not anymore possible such window is shifted vertically to underneath-left boundary position (e.g. after the first horizontal shifting movement the indices run over [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$m<y 2m$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$0<x n$$ as illustrated in Fig. (see window III). GLCM algorithm terminates once none of the above stjpg is further possible (i.e. W is too close to the right bottom corner of the I ). Every obtained W is the subject of further calculations. Co-occurrences of pixel values [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i=0, , q-1$$ are noted for a direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = 0^{ }, 45^{ }, 90^{ }, 135^{ }$$ and maximal distance d . Let [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{11}, ,w_{mn}$$ be coefficients of a given matrix W . Iterating through all of the coefficients one analyzes spatial correlations in their neighborhood. Currently analyzed pixel (coefficient of W in certain iteration) is named pixel of interest p . Definition of directions for maximal distance [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$d=2$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$p = w_{42}$$ is presented in Fig. . Assume p is [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{42}=0$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ =0^{ }$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$d=2$$ . Rows in [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ store the value of p and columns values of co-occurrence pixels being in relation to p . Since [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{42}=0$$ , the only coefficients in the first row of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ can be incremented. Then, based on Fig. , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{43}$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{44}$$ should be now analyzed. Coefficient [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{43}=0$$ implies that value in [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ in the first column ( [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}[1,1]$$ ) is incremented. Analogously, for [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{43}=3$$ [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}[1,4]$$ is increased. Both stjpg are presented in Fig. . In the same manner calculations for every W are performed analyzing all of the [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{ij}$$ for a direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ . Algorithm’s result is four matrices: [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{0^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{45^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{90^{ }}$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{135^{ }}$$ whose coefficients are then converted to the probabilities by dividing each of its values by mn . In the next step, the following statistical measures called GLCM features based on coefficients in matrices [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }$$ can be computed: Contrast, Correlation, Energy, Homogeneity, Autocorrelation, Cluster Prominence, Inverse Difference and Dissimilarity . Noticeably each of these features is calculable separately for each [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = 0^{ }, 45^{ }, 90^{ }, 135^{ }$$ . Having computed [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{0^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{45^{ }}$$ , [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{90^{ }}$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy^{135^{ }}$$ their mean [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Energy = 0.25(Energy^{0^{ }}+Energy^{45^{ }}+Energy^{90^{ }}+Energy^{135^{ }})$$ is next derived. Originally GLCM based features are calculated on the grayscale image i.e. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_t^{gray}$$ . A further improvement is obtained here once GLCM features are computed separately for each of the RGB image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_t^{RGB}$$ channel yielding in total 24 features. In contrast to the work by Duan et al. , which proposes an extension that utilizes images in a different color space than RGB for GLCM feature computation, the experiments incorporating different color spaces for GLCM feature computation for microscopic images of microorganisms do not result in an improvement of classification accuracy. Indeed, the best results presented in Table are reached for standard RGB color space, whereas other color channels might decrease classification accuracy from 87.4% to 74.6% using Y’IQ color space and 50 times repeated 10-Fold cross-validation technique and ELM with mish activation function and 1000 neurons on Dataset1.   In general GLCM is not the only method of texture features extraction. Commonly GLRLM is used together with GLCM. GLRLM creates features that describe spatial relationships between pixels having values of a given length of co-occurrence sequences and follows the concept of a run-length encoding outlined below. Non-overlapping sliding of window W similarly to GLCM is performed on the image quantized to the q levels. GLRLM is defined as four matrices [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{ }$$ each of size [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q r$$ , where r is a maximal run-length of co-occurrence pixels. Rows of that matrix are described as gray level values [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i = 0, ,q-1$$ and columns as run-lengths [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r_j=1, ,r$$ . For every obtained W run-length of co-occurrence pixel values q for direction [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and maximal run-length r is computed. Assume [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{11}, ,w_{mn}$$ form the coefficients of particular W . The directions are similarly defined as in GLCM method (see Fig. ). In GLRLM the whole window W for each of the gray levels [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i$$ is analyzed. To illustrate the latter for a given W with [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = 45^{ }$$ a maximal run-length [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r = 3$$ and take [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i = 0$$ were set. Then number of run-length sequences of length [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r_j=1, r$$ is noted. Based on such choice, the values in [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{45^{ }}$$ are incremented e.g. there are 2 sequences of length 2 (means that some row in 2nd column should be enlarged) for gray level equal 0 (first row), so coefficient of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^{45^{ }}[1,2]$$ is increased by 2 (see Fig. ). Similarly, all other entries of matrix [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^$$ are incremented for every window W based on input image and for every gray level [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i$$ . Consequently all coefficients of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^$$ are divided by qr . As a result matrix [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$M^$$ contains now probabilities of co-occurrence gray levels [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$q_i$$ for run-length [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$r_j$$ and maximum run-length r for a given angle [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ . Based on computed probabilities the following statistical measures for GLRLM (see Zhou et al. ) are calculated: Short Run Emphasis (SRE), Long Run Emphasis, Grey Level Non-uniformity, Run Length Non-uniformity, Run Percentage, Low Grey Level Run Emphasis, High Grey Level Run Emphasis, Short Run Low Gray Level Emphasis, Short Run High Gray Level Emphasis, Long Run Low Gray Level Emphasis and Long Run High Gray Level Emphasis . Similarly to GLCM, each of those features is calculated for four different directions [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and the final value of them is defined as e.g. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$SRE = 0.25(SRE^{0^{ }}+SRE^{45^{ }}+SRE^{90^{ }}+SRE^{135^{ }})$$ . Experiments to improve classification accuracy based on GLRLM features are conducted not to convert image to grayscale, but to calculate GLRLM on each of the color channels separately. Unfortunately, the same approach in case of GLRLM did not improve the results. A comparison is performed between the use of GLRLM feature generation based on images [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I_c$$ . Surprisingly, better results are achieved using image [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$I^{RGB}_t$$ , whereas conducting the same experiment on GLCM do not lead to further improvement in the results.   In total for each dataset the same set of features organized into three color spaces: RGB, I1I2I3, and XYZ (as explained in details in section before). The following abbreviations apply: RGB: Red (R), Green (G), Blue (B) I1I2I3: Intensity channels I1, I2, I3 XYZ: X, Y, Z color space Each color channel has the following statistical features: I1, I2, I3: Mean, Std, Mean (H), Std (H), Kurtosis (H), Skewness (H), Entropy (H), Energy (H) X, Y, Z: Mean, Std, Kurtosis, Skewness GLCM (Gray-Level Co-occurrence Matrix) features: Contrast, Correlation, Energy, Entropy, Homogeneity, Autocorrelation, Cluster Prominence, Inverse Difference, Dissimilarity (for R, G, B channels) GLRLM (Gray-Level Run Length Matrix) features: SRE, LRE, GLN, RLN, RP, LGRE, HGRE, SRLGE, SRHGE, LRLGE, LRHGE Notably the mathematical equations behind these metrics (e.g., Mean, Kurtosis, Skewness, GLCM, GLRLM) can be obtained from corresponding well-known machine vision literature , .   The ELM network is a fully connected feed-forward neural network (see Huang et. al. ) designed to solve classification problems. This network consists of one hidden-layer with McCulloch-Pitts neurons . The precise optimal number of hidden-layer neurons is each time evaluated empirically as there is no up-front theoretical method determining such value to maximize the classification accuracy. ELM is characterized by a simple topology structure. In network learning, the utilization of the Moore-Penrose pseudoinverse operation generally involves the application of singular value decomposition (SVD). The latter permits to shorten a time-consuming iterative process of adjusting values of network’s weights. Assume that after image segmentation, feature calculation and feature selection dataset consists of N images. It can be described as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = \{(x_i,t_i)\}_{i=1}^{N}$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$x_i$$ is i -th vector of d features and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ is i -th label of class to which selected image belongs. Then input layer of ELM has d neurons, defining also input matrix X introduced in . The number of output layer neurons coincides with the number of image classes M determining the classification task. Thus the output layer is numerically represented by a vector [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y = \{y_i\}_{i=1}^{M}$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_i  {R}$$ . Note that here the maximal value in vector y at index [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_{max}$$ is understood as assigning an input image to the [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_{max}$$ class. In order to test the ELM network, the i -th class label [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ needs to be appropriately formatted and compared to y . A possible approach is to format [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$\{t_{ij}\}_{j=0}^m$$ that has all values set to zero except one element at p -th index that is set to one. Then [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i$$ can be written as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i = [0, ,1, ,0]$$ where p -th element indicates that i -th image affiliates to p -th class of images ( [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_{i_p}=1$$ ). Only the greatest value of [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_i$$ at index [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$y_i^{max}$$ is chosen. Correct classification is observed if and only if [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$t_i[y_i^{max}] = 1$$ for i -th picture. This concept of formatting is called One-Hot Encoding or [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$1-of-K$$ scheme . Let now L be a number of units in single hidden layer network that is chosen empirically. Since ELM is a fully connected network, weights between input layer and hidden one are represented by a matrix W from , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$w_{ij}$$ denotes the weight between i -th neuron in input layer and j -th neuron of hidden layer. The corresponding bias values are assumed to be represented by a vector [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$b = \{b_i\}_{i=1}^{N}$$ . In ELM network both W and b are filled with values that are randomly chosen from the interval [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$[-1,1]$$ using a uniform distribution. Then the outputs of hidden units are calculated as expressed in matrix H (see , where f is any continuous function [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$f:  {R}  {R}$$ . 1 [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ X =  x_{11} &  & x_{1N}\\  &  &  \\ x_{d1} &  & x_{dN} ,  W =  w_{11} &  & w_{1L}\\  &  &  \\ w_{d1} &  & w_{dL} ,  H =  f( _{i=1}^{d}x_{i1}w_{i1} + b_1) &  & f( _{i=1}^{d}x_{i1}w_{iL} + b_1)\\  &  &  \\ f( _{i=1}^{d}x_{iN}w_{i1}+ b_N) &  & f( _{i=1}^{d}x_{iN}w_{iL} + b_N) . $$ Weights [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ between hidden and output layers are still unknown. In order to obtain them, the equitation [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ has to be evaluated. Indeed since H with probability equal to 1 is irreversible and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$\| H -Y \| = 0$$ (see Huang et al. ), the system [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ cannot be directly solved. Instead, one finds the optimal parameter [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ = arg\ \| H -T \| ^2 = H^ T$$ , where [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ is the Moore-Penrose generalized inverse of H (see Rao and Mitra ) and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$T = {t_i}_{i=1}^{N}$$ describing image appurtenance to the given class formed using concept of one hot encoding mentioned in the above paragraph. Moore-Penrose pseudoinverse of a real matrix H must satisfy following conditions: 1. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$HH^ H = H,$$ 2. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^ HH^ = H^ ,$$ 3. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$(HH^ H)^T = HH^ ,$$ 4. [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$(H^ H)^T = H^ H.$$ Pseudoinverse is known to be unique and pseudoinverse of a non-singular quadratic matrix is identical to the ordinary inverse. For matrices that are singular pseudoinverse solves [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y = H$$ in terms of least-squares problem. Calculating [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ gives solution [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ so that [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H$$ is as close as possible to Y in terms of Euclidean distance. Note here that H generically corresponds to either overdetermined or underdetermined system. Additionally, there can be many methods to calculate [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ , one of them recalled here is orthogonal projection method (see Rao and Mitra ). If [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$L>d$$ (number of hidden units is greater than number of input units) then system [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ is underdetermined and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^ = H^T(HH^T)^{-1}$$ , otherwise system [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$Y=H$$ is overdetermined and [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^ = (H^TH)^{-1}H^T$$ . In practical applications of ELM [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$L d$$ eventuates. In classification tasks L usually totals a few thousands. Then, operation of inverse [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$HH^T$$ could be time-consuming especially taking into account the fact that operation of matrix pseudoinverse cannot be parallelized. To improve learning speed of network QR matrix decomposition is used (see Strang ). QR decomposition of matrix A is defined as [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$A=QR$$ , where R is upper triangular and Q is orthogonal ( [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$QQ^T =  {I}$$ ). Using QR decomposition, the matrix [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$H^$$ can be computed according to: [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$ H^{ } = H^T(HH^T)^{-1} = R^TQ^T(R^TQ^TQR)^{-1} = R^TQ^T(R^TR)^{-1} = Q^TR^{-1}. $$ The ELM network gives a solution that minimizes mean square error between predicted [12pt]{minimal} 				 				  				  				  				 				 				 				 				$$$$ and actual Y . Thanks to assigning random values to weights and bias between input and hidden layer, the trained network is very resistant to overtraining. The computed solution is not a local minimum of optimization criterion like in MLP with backpropagation algorithm based on gradient descent , . The network training process is much faster than other Machine Learning methods like MLP, RF and SVM. The main drawback of this process relies on finding L and f that maximize classification accuracy.   