CA23125 - The mETamaterial foRmalism approach to recognize cAncer (TETRA); Working Groups 2, 2026 Mathematical modelling of highly disordered anisotropic structures. Part 2.3. Classification algorithm for automatic detection for cancerous zones in tissue. Vladimir Mityushev July 23, 2026 Abstract The introduction to selected aspects of DIP (Digital Image Processing) is presented. Further we employ this knowledge to construct image process- ing pipelines for specific images. We illustrate all concepts using Python programming language. As a short introduction, we can recommend an ap- pendix of our book [17]. You can also study free Scientific Python Lectures [3]. Complete examples of the code from this part can be downloaded from [1]. You can run it and experiment. We use the OpenCV (Open Computer Vision) library for image processing. We learn the specific algorithms in the library that are useful for image processing of tissues and selecting parts suspected to be cancer cells. We focus on cells in medical applications. The developed aRVE methodological, see Parts 2.1 and 2.2, and practical applica- tion enables a systematic treatment of the variability in random composites within the framework of homogenization principles and asymptotic analysis. Application to the diagnosis and classification of glioma structures requires a set of corresponding pictures to apply the developed machine learning classi- fication algorithm. To summarize, the theoretical basis and algorithms have been developed. In order to go further to practical problems we need: • physical constants of the image components; • a number of healthy and diseased samples to process their and use for the classification of new pictures by supervised machine learning. 1
1 Digital Image Processing (DIP) We begin with the main concepts of DIP concerning thresholding, detection of zones and their classification. They will be used in practical applications. We want to be manageable due to a specific aim in our head, and since DIP is a vast discipline, there are already excellent resources that treat this topic in great generality. You can consult [12] for very precise mathematical exposition or more concise yet still extremely precise [16]. We will illustrate concepts constructing Python code with the help of OpenCV library [2], a robust and well-optimized set of image processing tools for real-world applications. The library is successively amended by new algorithms and methods and is a de facto standard in Image Processing. The library is written in C language due to optimization and possible applications in embedded systems. However, a version for Python language allows us to use the benefits of high-level language if we do not need real-time responses. There are many excellent and practical books that introduce to the library, e.g., [21, 13], and reference books, e.g., [15]. The package Mathematica is also used in DIP. It is worth noting that Python and Mathematica complement each other to construct algorithms for automatic detection for zones in tissue. In the following subsections, you will be introduced to all necessary topics from DIP to gain working knowledge for processing issues images, materials images, and biological microscopic images. 1.1 Thresholding Having a grayscale image, we want to extract some shapes of a given inten- sity level. This can be done by selecting pixels above the fixed (intensity) threshold. We can use histograms to specify the threshold or use automatic methods. The general idea of thresholding θt is to modify pixel-wise the value de- pending on the value of threshold t and maximal value V (typically the maximal intensity V = 255). In a simple binary thresholding we have θV t (f )(x, y) = 0 if f (x, y) < t V if f (x, y) ≥ t (1.1) We can also invert the values. This is the most useful in creating masks. First, we create an image that contains two types of granules/circles of different sizes. We, moreover, add a disturbance to the background. We then convert the image to the grayscale. In real life, the complexity is larger 2
1 import numpy as np 2 import matplotlib.pylab as plt 3 import cv2 4 5 #prepare figure 6 image = np.zeros((512,512,3), np.uint8) 7 image = image + 255 8 9 #perturbation 10 cv2.rectangle(image,(100,0),(500,400),(240,240,240),-1) 11 cv2.rectangle(image,(0,100),(400,500),(0,240,240),-1) 12 13 #generate class 1 14 x = 512* np.random.rand(200) 15 y = 512* np.random.rand(200) 16 for i,j in zip(x,y): 17 cv2.circle(image,(int(i),int(j)), 5, (100,100,100), -1) 18 19 #generate class 2 20 x = 512* np.random.rand(50) 21 y = 512* np.random.rand(50) 22 for i,j in zip(x,y): 23 cv2.circle(image,(int(i),int(j)), 20, (10,10,10), -1) 24 25 imageGray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) 26 cv2.imwrite("output/initialImageThresholding.png", imageGray) which produces Figure 1. The next is to construct histogram of intensity distribution 1 histGray= cv2.calcHist([imageGray],[0],None,[256],[0,256]) 2 plt.plot(histGray/max(histGray)) 3 plt.xlabel("bin") 4 plt.ylabel("normalized cout") 5 plt.savefig("output/GrayHistogramThresholding.png") 6 plt.show() that gives Figure 2. We can note in Figure 2 the bin close to zero repre- senting bigger darker balls in Figure 1, and the second bin closer to 100 that represents smaller, lighter balls. These two kinds of balls are granules we 3
Figure 1: Initial image converted to grayscale. Figure 2: Intensity histogram for image 1. 4
Figure 3: Extraction of larger granules from Figure 1. want to extract. Moreover, we have three bigger peeks above 200: two light rectangles of big area and the last one at 255 represents the background white color. The first kind of granules can be extracted easily by setting the threshold value to, e.g., 50. The pixels equal to or above this value will be set to white. We can do it using the following code 1 _ , mask1 = cv2.threshold(imageGray, 50, 255, cv2.THRESH_BINARY) 2 cv2.imwrite("output/MaskThresholding1.png", mask1) which produces Figure 3, which can be used as a mask for selecting the first kind of balls. For selecting the second kind of balls we have to subtract from Figure 1 the balls from image 3. That can be done by negating the mask 3 and using exclusive OR (XOR) binary operation with the original image 1. Then we can do thresholding below, e.g., 150 1 _ , mask2 = cv2.threshold(cv2.bitwise_xor(imageGray, 2 cv2.bitwise_not(mask1)), 150, 255, cv2.THRESH_BINARY) 3 cv2.imwrite("output/MaskThresholding2.png", mask2) that produces Figure 4. This figure can also be used as a mask for extracting smaller balls. One can observe that some of these balls were shadowed by 5
Figure 4: Extraction of smaller granules from Figure 1. the darker, bigger ones, and therefore, they are only partially restored. We can remove these parts using morphological operations as described above, however, we will not do it here. At this point we have desired segmentation of the image into two classes of granules. The same extraction can be done using cv2.inRange function that is supplied with the boundaries of the intensity interval from which we extract pixels as in the following code 1 mask2 = cv2.bitwise_not(cv2.inRange(imageGray, 50, 150)) 2 cv2.imwrite("output/MaskThresholding2.png", mask2) We also inverted colors to have granules black and color white (opposite to the colors used in a mask). The function cv2.inRange can also be used to threshold color images, and then the boundary points are BGR vectors for minimal and maximal values of each color channel. We can now focus on examining automatic thresholding. We first examine Otsu’s method [18], which tries to split pixels into different classes and tries to maximize in-between-class variance, see, e.g., [12]. It is realized, but the following code 1 _, imageOtsu = cv2.threshold(imageGray, 0, 255, cv2.THRESH_OTSU) 2 cv2.imwrite("output/MaskThresholdingOtsu.png", imageOtsu) 6
Figure 5: Otsu’s method applied to Figure 1. that produces Figure 5 that removes background. We can also use adaptive thresholding using the mean local value of pixels minus a constant. We have to provide the size of the local neighborhood. This is realized by 1 imageAdaptiveMean = cv2.adaptiveThreshold(imageGray,255, 2 cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY,51,1) 3 cv2.imwrite("output/MaskThresholdingAdaptiveMean.png", 4 imageAdaptiveMean) that produces Figure 6. The last two parameters in cv2.adaptiveThreshold are the size of the local neighborhood and the constant C that is subtracted from the local mean, giving the local threshold. One can note that some artifacts resulting from removing background rectangles are present. We can also apply Gaussian adaptive thresholding, which determines the threshold as Gaussian-weighted local neighborhood values minus a constant. This algorithm is applied in the following code 1 imageAdaptiveGaussian = cv2.adaptiveThreshold(imageGray,255, 2 cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY,51,1) 3 cv2.imwrite("output/MaskThresholdingAdaptiveGaussian.png", 7
Figure 6: Adaptive mean thresholding applied to Figure 1. 4 imageAdaptiveGaussian) that produces Figure 7. In this case, also, some artifacts are visible. Summing up, the basic automatic methods are less helpful in extracting granules in such a simple example. 1.2 Contours Having extracted granules, we want to characterize them. In this section we will learn how to extract contours, fit geometric shapes to them, and compute moments of contours. First, we must prepare the image. We use an academic example of random balls. The difficulty will be with overlapping balls. The following code 1 import numpy as np 2 import pandas as pd 3 import matplotlib.pylab as plt 4 import cv2 5 6 #prepare figure 7 image = np.zeros((512,512,3), np.uint8) 8 image = image + 255 8
Figure 7: Adaptive gaussian thresholding applied to Figure 1. 9 10 #generate granules 11 np.random.seed(1) 12 x = 512* np.random.rand(50) 13 y = 512* np.random.rand(50) 14 for i,j in zip(x,y): 15 cv2.circle(image,(int(i),int(j)), 20, (255,0,0), -1) 16 17 cv2.imwrite("output/initialImageContours.png", image) creates Figure 8. First, we create a grayscale mask 1 imageGray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) 2 imageGray = cv2.bitwise_not(imageGray) 3 cv2.imwrite("output/maskImageContours.png", imageGray) that gives Figure 9. Since some granules are merged we can try to disconnect them using erosion multiple times. However, overlapping is too big in some cases to separate them by erosion without significantly removing white pixels from well-separated images. Moreover, in real applications, granules have 9
Figure 8: Initial image for detecting contours. Figure 9: Grayscale image of Figure 8. 10
Figure 10: Detected contours (red) from Figure 8. irregular shapes and are connected to such an extent that they are impossible to distinguish even by the human eye. OpenCV has many specialized methods that can detect contours. Con- tours are lists of points that sit at the edge of shapes. In the simplest case, the algorithm tries to find the contours(borders) of connected components. We use a general one for detecting contours on grayscale image 1 contours, _ = cv2.findContours(imageGray, cv2.RETR_LIST, 2 cv2.CHAIN_APPROX_TC89_L1 ) 3 #draw contours 4 imageGrayContours = image.copy() 5 #imageGrayContours = np.zeros(image.shape)+255 6 cv2.drawContours(imageGrayContours, contours,-1,(0,0,255),3) 7 cv2.imwrite('output/contoursImageContours.png', 8 imageGrayContours) which gives Figure 10. One can observe that some blobs are treated as a single connected component. The function cv2.findContours has various options and types of re- turned data. For detecting granules we need only a flat structure of the list (option cv2.RETR LIST). In general, it can return a hierarchical struc- ture of contours using cv2.RETR TREE instead. We can also use various ap- 11
proximation methods that optimize the number of points. Since we need many points in order to match ellipses to the contours, we used option cv2.CHAIN APPROX TC89 L1 which uses the Teh-Chin chain approximation algorithm [20]. We can use other Teh-Chin algorithm by option cv2.CHAIN APPROX TC89 KCOS, or simple approximation with a small number of points cv2.CHAIN APPROX SIMPLE, or no approximation cv2.CHAIN APPROX NONE - stores all points of the con- tours. The contour is drawn on the image imageGrayContours all (−1 contours of all indices) detected contours using a red line of thickness 3. Since we know that the contours are balls, we can apply more specialized algorithms to detect circles that are available in OpenCV. We use the method based on Hough transform [8] 1 img = imageGray.copy() 2 colorimg = image.copy() 3 circles = cv2.HoughCircles(img,cv2.HOUGH_GRADIENT, 1, 20, 4 param1=20, param2=10, minRadius=2, maxRadius=50) 5 circles = np.uint16(np.around(circles)) 6 for circle in circles[0,:]: 7 # draw the outer circle 8 cv2.circle(colorimg,(circle[0],circle[1]),circle[2], 9 (0,0,255),2) 10 # draw the center of the circle 11 cv2.circle(colorimg,(circle[0],circle[1]),3,(0,0,0),3) 12 cv2.imwrite('output/HoughCirclesDetectorImageContours.png', 13 colorimg) that produces Figure 11. One can note that more balls are detected, however, the algorithm does not spot some overlapping granules. The list of circles returned by the cv2.HoughCircles method represented by x and y coordinates (indices 0 and 1) of the center and the radius (index 2). Let us return to contours, as they are more general types of closed curves used for bonding regions in the images. We want to draw around the con- tour bounding rectangles. Moreover, we cut the rectangular image enclosing the contours and save them in images, moreover, we enumerate detected contours. This is realized by the following code 1 img = image.copy() 2 for idx, contour in enumerate(contours): 3 12
Figure 11: Hough algorithm for detecting circles applied to Figure 8. 4 #draw minimum area rectangle 5 minRect = cv2.minAreaRect(contour) 6 box = cv2.boxPoints(minRect) 7 box = np.int0(box) 8 cv2.drawContours(img, [box], -1, (0,0,255), 1) 9 #enumarate contours 10 cv2.putText(img, '{}'.format(idx), (box[0][0], box[0][1]), 11 cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0,0,0), 1, cv2.LINE_AA) 12 13 #extract image from bounding rectangle and write to files 14 x, y, w, h = cv2.boundingRect(contour) 15 cropped = image[y:y + h, x:x + w, :].copy() 16 feature = cv2.resize(cropped, (20, 20)) 17 cv2.imwrite("./output/figs/{}.png".format(idx), feature) 18 cv2.imwrite('output/MinAreaRectangleImageContours.png', img) which produces Figure 12. Extracted images can be used for (unsupervised) training Machine Learning algorithms to discriminate their types if there are more types of them. We will return to Machine Learning later in this part. When the granules have no circular symmetry, we sometimes want to know their orientation and non-symmetry. The easiest way is to fit ellipses and extract their parameters. OpenCV has a function that can fit an ellipse 13
Figure 12: Minimal area rectangle for detected contours from Figure 10. using at least five points of contours using the least square method. We do it in the following code 1 img = image.copy() 2 ellipses = [] #list for storing ellipsis parameters 3 for idx, contour in enumerate(contours): 4 if contours[idx].size < 10 : 5 continue 6 ellipse = cv2.fitEllipse(contour) # ellypse is ((x,y), 7 (majorAxis, minorAxis), angle) 8 #remove bigger ellipses 9 if max(ellipse[1]) < 200: 10 ellipses.append({"id": idx, "x":ellipse[0][0], 11 "y":ellipse[0][1], "majorAxis": ellipse[1][0], 12 "minorAxis": ellipse[1][1], "angle": ellipse[2]}) 13 cv2.ellipse(img, ellipse, (0,0,255), 3, cv2.LINE_AA) 14 cv2.putText(img, '{}'.format(idx), (int(ellipse[0][0]), 15 int(ellipse[0][1]) ), cv2.FONT_HERSHEY_SIMPLEX, 16 0.2, (0, 0, 0), 1, cv2.LINE_AA) 17 #save parameters of ellipses 18 df_ellipses = pd.DataFrame.from_records(ellipses) 19 df_ellipses.to_csv("output/ellipses.csv") 14
Figure 13: Fitted ellipses to contours from Figure 10. 20 #write image with elipses 21 cv2.imwrite('output/EllipsesImageContours.png', img) that produces the image in Figure 13. The ellipse is characterized by the center (x, y), major axis, minor axis, and orientation given by the angle between the positive 0X axis and the major axis. In the code, these data are packed into a Python DataFrame structure and saved in ellipses.csv file1 Sometimes, we want to characterize contours more precisely. We can use moments of distribution. For a contour c(x, y) we can get static real moments (sometimes called ’spatial or raw moments’) mij = X x,y c(x, y)xiyj , (1.2) central moments μij = X x,y c(x, y)(x − ¯x)i(y − ¯y)j , (1.3) 1CVS (comma-separated values) is a simple and portable file format for saving tabular data in a text file. Each row of the table corresponds to the row of the file, and columns are represented by data separated by a comma or other separator. These files can be easily imported into a spreadsheet program. 15
where the ’mass center’ is characterized by ¯x = m10 m00 , ¯y = m01 m00 , (1.4) and normalized central moments (ν =′′ nu′′) defined as νij = μij m(i+j)/2+1 00 . (1.5) Since μ00 = m00, ν00 = 1, ν10 = μ10 = ν01 = μ01 = 0, so these values are not returned. The considered static real moments correspond to the static complex moments introduced in Part 2.2. The complex Green’s formulas allow us to treat line integral over a closed contour as a surface area over the area enclosed by the contour. The same holds for the static real moments. The contours with self-intersections can have zero m00 even if the area is not zero due to the implicit orientation of areas in Green’s formula. Moreover, the value depends on rasterization, and therefore, the computed values can be treated only as approximations. Moreover, the more points of the contour yield a more accurate approximation. The contours (viewed as continuous curves/before discretization) can be transformed by scaling, translations, and rotations (after undoing discretiza- tion on the raster lattice). We are interested in having the same moments for two of the same contours, even if one is a transformed version of the second one. Therefore, in analyses, we should use moments that are invariant under specific symmetry groups. A well-known fact is that the central moments are invariant with respect to translations. The moments νij are invariant with respect to translations as defined by central moments, and also with respect to scaling. The natural next step is to define the invariant moments concerning scaling, rotations, and translations. These seven moments were defined by Hu [14], and then amended [9] and generalized for spatial objects 16
[10]. Hu moment invariants are defined as I1 = ν20 + ν02, I2 = (ν20 − ν02)2 + 4ν2 11, I3 = (ν30 − 3ν12)2 + (3ν21 − ν03)2, I4 = (ν30 + ν12)2 + (ν21 − ν03)2, I5 = (ν30 − 3ν12)(ν30 + ν12)[(ν30 + ν12)2 − 3(ν21 + ν30)2]+ +(3ν12 − ν03)(ν21 + ν03)[3(ν30 + ν12)2 − (ν21 + ν03)2], I6 = (ν20 − ν02)[(ν30 + ν12)2 − (ν21 + ν03)2] + 4ν11(ν30 + ν12)(ν21 + ν03), I7 = (3ν21 − ν03)(ν30 + ν12)[(ν30 + ν12)2 − 3(ν12 + ν03)2]− (ν30 − 3ν12)(ν21 + ν03)[3(ν30 + ν12)2 − (ν21 + ν03)2]. (1.6) They are not all independent. Moreover, I7 is reflection antisymmetric and the remaining ones reflection symmetric, they can be used to detect mirror symmetry of the contour. In addition, there are more invariants of fourth order [9]. The Hu moment invariants can be computed in OpenCV using cv2.HuMoments function using the result of cv2.moments that computes mij , μij and νij moments. The following code, being a continuation of the previous code, computes moments and Hu moment invariants and saves it to a CSV file 1 import math 2 huMomentsList = [] 3 for idx, contour in enumerate(contours): 4 #calculate moments 5 moments = cv2.moments(contour) 6 print("contour nr ", idx, ", moments:", moments) 7 #calculate Hu moments 8 huMoments = cv2.HuMoments(moments) 9 #scaling Hu moments 10 for i in range(0,7): 11 if huMoments[i] == 0: 12 huMoments[i] = 0.0 13 else: 14 huMoments[i] = -1* math.copysign(1.0, huMoments[i]) 15 * math.log10(abs(huMoments[i])) 17
Figure 14: Hu moments for contours from Figure 10. 16 17 huMomentsList.append({'h0':huMoments[0], 'h1':huMoments[1], 18 'h2':huMoments[2], 'h3':huMoments[3], 19 'h4':huMoments[4], 'h5':huMoments[5], 20 'h6':huMoments[6]}) 21 22 df_HuMoments = pd.DataFrame.from_records(huMomentsList) 23 df_HuMoments.to_csv("output/HuMoments.csv") Since Hu moment invariants can have arbitrary values, we adjust their range using log10 transform and preserve their sign. It is standard operation, and it can be helpful when we use moments for Mchine Learning algorithms that require data to be properly normalized/standardized to have a similar range. This code generates the file presented in Figure 14. In the present section, the discussed static moments are used for the image analysis. It is worth noting that the static moments affect the physical properties in combination with structural sums. This is explicitly shown by 18
the formulas from Part 2.1 and 2.2. 1.3 Classifying Once we ’segmented’ image into granules and identify potentially interesting ones, the next step is to identify them. The most reliable is the ’hand-made’ classification by a skilled human expert. On the other hand, this is the slowest and least automatic way of classification. In some important cases where the error should be minimized, like cancer cell detection, the human expert should have the last word in such classification. However, automatic methods can help in such classification and speed up the process. In this part, we make a short review of Machine Learning (ML) methods that can be applied to the classification of images. We focus mostly on the methods where no expert insight is available (unsupervised learning). We also explain how to use other types of Machine Learning where expert data is present for learning before classification (supervised learning). The basics of ML in the context of mathematical modeling can be learned in our book [17]. Moreover, a more formal yet practical introduction to Machine Learning are the books [19, 11]. This part cannot be treated as an introduction to the vast area of ML. However, you will learn the idea and see an example that can be generalized to any other, more advanced ML method. Since computers are good only at crunching numbers, all data, before processing, must be converted into numbers or vectors of numbers. The general idea of ML methods is to distinguish separate classes in the data and associate class labels to them. We distinguish two main cases: • supervised learning - when class labels are present in the training data; Then ML algorithm learns on examples to associate labels to the new data on training. This type of algorithms focus on making precise borders between areas with points associated with different class labels (decision regions). In this class we have algorithms like SVM (Sup- port Vector Machines) and various architectures of Artificial Neural Networks. • unsupervised learning - when that data does not contain class labels, the ML algorithm must infer how many classes there are and label them. The classes are typically determined as clusters in the data space. All kinds of clustering algorithms like k-Means or Gaussian Mixture Model represent this class. Supervised learning can be made with the data that are labeled. If no such data are available for our case, the way to detect differences in data is 19
to use unsupervised learning. In the later part we will focus on unsupervised learning. Already trained supervised learning algorithms can be treated as a replacement in the code encountered in this part. Many ML algorithms are implemented in the Scikit Learn library [4]. The work with classifiers is rather standard and consists of two steps: 1. fit(...) - learn classifier; 2. predict(...) - use a classifier to predict classes; usually used for new data not used in the learning process; A prominent example of unsupervised ML algorithms is the k-means method. It tries associating data to k clusters, minimizing WSSE (clus- ters within the sum of square error). The input is the number of clusters, however, it can be determined. In order to understand the idea, we present the following toy example, where we create artificial data with some clusters and then apply the k-Means algorithm. More details can be found in our book [17]. 1 import matplotlib.pylab as plt 2 3 #preparing data 4 from sklearn.datasets import make_blobs 5 X, y = make_blobs(n_samples=1000, centers=3, n_features=3, 6 random_state=0) 7 8 # plot data 9 plt.scatter(X[:,0], X[:,1], c=y) 10 plt.xlabel("$x1$") 11 plt.ylabel("$x2$") 12 plt.savefig("output/data-kMeans.png") which produces three clusters in Figure 15. Three clusters are clearly visible, and we expect to distinguish them. The easiest way to determine the number of clusters is to plot WSSE vs the number of clusters (the scree plot) and identify characteristic elbow. The following code 1 from sklearn.cluster import KMeans 2 wsse = [] 3 for k in range(1, 20): 4 kmeans = KMeans(n_clusters=k, n_init=1) 20
Figure 15: Initial data for k-Means algorithm. 5 kmeans.fit(X) 6 wsse.append(kmeans.inertia_) 7 8 #make plot 9 plt.plot(range(1, 20), wsse) 10 plt.xlabel("Clusters") 11 plt.ylabel("WSSE") 12 plt.savefig("output/kMenas_Scree.png") gives Figure 16. One can see ’elbow’ with ’joint’ at k = 3 clusters. Above this number of clusters, increase does not improve WSSE much. We can now make k-Means clustering using k = 3 clusters 1 kmeans = KMeans(n_clusters=3, n_init=1) 2 kmeans.fit(X) 3 ##clusters labels 4 labels = kmeans.predict(X) 5 ##locations of centroids 6 centroids = kmeans.cluster_centers_ 7 8 #plot 9 plt.scatter(X[:,0], X[:,1], c=labels) 10 plt.scatter(centroids[:,0],centroids[:,1], c='r', marker='v', 11 label="centers") 12 plt.legend() 13 plt.xlabel("x1") 21
Figure 16: Scree plot. 14 plt.ylabel("x2") 15 plt.savefig("output/kMeans-clustering.png") 16 plt.show() that gives plot in Figure 17. The above example was an ideal one. Real-world data are typically not so well separable. The next example will present the k-means algorithm for image classifi- cation, which will be useful in the next part of this part. We use handwritten digits [5] as input data. They are loaded by the following code 1 import numpy as np 2 import matplotlib.pylab as plt 3 4 from sklearn.datasets import load_digits 5 6 digits = load_digits() 7 X = digits.data 8 Y = digits.target 9 10 #reverse color 11 X = 255-X 12 13 #example digit 14 plt.imshow(X[100].reshape(8,8)) 15 plt.savefig("output/digitSingleOne.png") 16 plt.show() 22
Figure 17: Clustering result. Figure 18: An example of handwritten digit. that produces Figure 18. The image is encoded as a vector of 64 entries and has to be converted to 8 × 8 grayscale array. The next step is to estimate the number of clusters 1 from sklearn.cluster import KMeans 2 k = range(1,20) 3 WSSE = [] 4 for i in k: 5 model = KMeans(n_clusters = i, n_init='auto') 6 model.fit_predict(X) 7 WSSE.append(model.inertia_) 8 23
Figure 19: Scree plot for handwritten digits data. 9 plt.plot(k,WSSE) 10 plt.xlabel('K') 11 plt.ylabel('Sum of Squared Errors') 12 plt.savefig("output/digitsScree.png") 13 plt.show() 14 that gives Figure 19. One can note that the elbow joint is about 10, which is expected. Therefore in the next step we can use the k-Means classifier to cluster digits into similar sets 1 clusters = 10 2 kmeans = KMeans(n_clusters=clusters, n_init='auto') 3 kmeans.fit(X) 4 Z = kmeans.predict(X) 5 6 #print selected class 7 ##get images of numbers from given class 8 clusterIndex = 2 9 #get images with a specific number of cluster 10 imageIndex = np.where(Z==clusterIndex)[0] 11 #plt.figure(figsize=(10,10)) 12 for k in range(0, imageIndex.shape[0]): 13 image = X[imageIndex[k],:] 14 image = image.reshape(8, 8) 15 plt.imshow(image, cmap='gray') 16 plt.savefig("output/digits/{}_{}.png". 24
17 format(clusterIndex, k)) plt.show() The results assigned to the class 2 can be inspected visually. It is not perfect, however, many images represents the same digit. Instead of the k-means algorithm, one can use any other unsupervised ML algorithm or trained, supervised ML algorithm. 1.4 Segmentation using Means algorithm One of the standard techniques of segmentation of an image is related to similar colors of objects we want to extract. If this is the case, we can do clusterization in color space, and then extract objects belonging to the same cluster (similar color). Clusterization can be performed using an arbitrary unsupervised ML algorithm. We present an example that uses the k-means algorithm. We start by preparing an image 1 import numpy as np 2 import pandas as pd 3 import matplotlib.pylab as plt 4 import cv2 5 6 #prepare figure 7 image = np.zeros((512,512,3), np.uint8) 8 image = image + 255 9 10 #perturbation 11 cv2.rectangle(image,(100,0),(500,400),(240,240,240),-1) 12 cv2.rectangle(image,(0,100),(400,500),(0,240,240),-1) 13 14 #set reproducibility of initial image 15 np.random.seed(1) 16 17 #generate class 1 18 x = 512* np.random.rand(200) 19 y = 512* np.random.rand(200) 20 for i,j in zip(x,y): 21 cv2.circle(image,(int(i),int(j)), 5, (100,100,100), -1) 22 23 #generate class 2 24 x = 512* np.random.rand(50) 25
Figure 20: Initial image. 25 y = 512* np.random.rand(50) 26 for i,j in zip(x,y): 27 cv2.circle(image,(int(i),int(j)), 20, (10,10,10), -1) that produces an image from Figure 20. Next, we localize the optimal number of classes 1 from sklearn.cluster import KMeans 2 3 k = range(1,5) 4 WSSE = [] 5 6 for i in k: 7 model = KMeans(n_clusters = i, n_init='auto') 8 model.fit_predict(df) 9 WSSE.append(model.inertia_) 10 11 plt.plot(k,WSSE) 12 plt.xlabel('K') 13 plt.ylabel('WSSE') 14 plt.grid(True) 15 plt.savefig("output/kmeansSegmentationColor.png") 26
Figure 21: Scree plot for Figure 20. 16 plt.show() that produces Figure 21. One can observe that 2 clusters are enough. The following code clusters data with respect to 2 clusters 1 NClusters = 2 2 3 km = KMeans( n_clusters= NClusters, n_init='auto') 4 labels = km.fit_predict(df) 5 #create mask from class labels 6 mask = labels.reshape(xShape, yShape, 1) 7 8 imgLabelled = [] 9 masks = [] 10 for label in range(NClusters): 11 maskClass = np.where(mask == label, 255, 0) 12 maskClass = maskClass.astype(np.uint8) 13 masks.append(maskClass) 14 cv2.imwrite("output/SegmentationColorMask{}.png". 15 format(label), maskClass) 16 imgLabel = cv2.bitwise_and(imageGray, imageGray, 17 mask = maskClass) imgLabelled.append(imgLabel) 18 cv2.imwrite("output/SegmentationColorLabel{}.png". 19 format(label), imgLabel) that produces segments in Figure 22, which were made using masks Figure 23 applied to Figure 20. We can see that the masks are inverse to each 27
Figure 22: Segmentated image for Figure 20. Figure 23: Masks that applied to Figure 20 result in Figure 22. 28
Figure 24: Masks for k = 4 clusters. other. It is instructive to consider a higher number of clusters. Figure 24 rep- resents masks for NClusters = 4 clusters in the k-Means algorithm. We observe that the k-Means algorithm separates small and big granules, as well as elements in the shaded region and its complement. 1.5 Summary Until now, we have presented many essential methods of image manipulation with the OpenCV library that can be used for extracting and classifying granules in an image. They mustn’t be used as a black box. They must be used in a specific order, and parameters must be adjusted to a particular situation/type of image. Some experience is needed, however, thanks to them you can automatize granulometry. 29
Figure 25: Initial image. Microscope image of cancerous mouse brain biopsy. 2 Examples In this section, we will describe more advanced examples that use methods from the previous part. We present four examples. 2.1 Example 1 In the first example we want to select darker granules from Figure 25. This is a microscope image of a tissue. In the first step we read the image and cut out the top part of the image 1 import numpy as np 2 import pandas as pd 3 import matplotlib.pylab as plt 4 import cv2 5 6 img1 = cv2.imread("Picture1.png") 7 8 #cut out printout 9 img1 = img1[50:-10,10:] 10 cv2.imwrite("output/Example1Cut.png", img1) 30
Figure 26: Cropped image that will be further processed. which gives Figure 26. We now improve contrast by histogram equalization 1 #histograms of colors 2 hist= cv2.calcHist([img1],[0],None,[256],[0,256]) 3 hist2=cv2.calcHist([img1],[1],None,[256],[0,256]) 4 hist3=cv2.calcHist([img1],[2],None,[256],[0,256]) 5 ##plot histograms 6 plt.subplot(1, 3, 1) 7 plt.plot(hist/max(hist)) 8 plt.title('Hist Blue : ') 9 plt.subplot(1, 3, 2) 10 plt.plot(hist2/max(hist2),color="green") 11 plt.title('Hist Green : ') 12 plt.subplot(1, 3, 3) 13 plt.plot(hist3/max(hist3),color="red") 14 plt.title('Hist Red : ') 15 plt.savefig("output/Example1ColorHistograms") 16 plt.show() 17 18 #histogram equalization 19 src = cv2.cvtColor(img1, cv2.COLOR_BGR2HSV) 20 h,s,v = cv2.split(src) 21 v = cv2.equalizeHist(v) 22 src = cv2.merge([h,s,v]) 31
Figure 27: Result of histogram equalization applied to Figure 26. 23 img1 = cv2.cvtColor(src, cv2.COLOR_HSV2BGR) 24 cv2.imwrite("output/Example1Equalization.png", img1) 25 26 #histograms of colors after normalization 27 hist= cv2.calcHist([img1],[0],None,[256],[0,256]) 28 hist2=cv2.calcHist([img1],[1],None,[256],[0,256]) 29 hist3=cv2.calcHist([img1],[2],None,[256],[0,256]) 30 plt.subplot(1, 3, 1) 31 plt.plot(hist/max(hist)) 32 plt.title('Hist Blue : ') 33 plt.subplot(1, 3, 2) 34 plt.plot(hist2/max(hist2),color="green") 35 plt.title('Hist Green : ') 36 plt.subplot(1, 3, 3) 37 plt.plot(hist3/max(hist3),color="red") 38 plt.title('Hist Red : ') 39 plt.savefig("output/Example1ColorHistogramsEqualization.png") 40 plt.show() that apart of color histograms, produces histogram-equalized Figure 27. In order to enhance contours, we apply bilateral filtering, convert to grayscale, and compute histogram in order to determine the threshold. The following code 1 #scale figure 32
Figure 28: Grayscale image after bilateral filtering to enhance gradients. 2 img1 = cv2.resize(img1,None,fx=4, fy=4, interpolation = 3 cv2.INTER_CUBIC) 4 # convert to BW 5 img1_BW = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY) 6 cv2.imwrite("output/Example1BW.png", img1_BW) 7 8 #bilateral filtering 9 img1_BW = cv2.bilateralFilter(img1_BW,9,75,75) 10 cv2.imwrite("output/Example1BWFiltering.png", img1_BW) 11 12 #compute gray histogram 13 hist1 = cv2.calcHist([img1_BW], [0], None, [256], [0,256]) 14 plt.plot(hist1.ravel()/hist1.max()) 15 plt.savefig("output/Example1HistogramBW.png") 16 plt.show() which produces Figs. 28, 29. One can note that the black dots are rep- resented by the pixels with intensity below 40. We use this information in thresholding and below. We also do morphological operations to delete small granules in the following code 1 ret1m, thresh1 = cv2.threshold(img1_BW, 40, 256, 2 cv2.THRESH_BINARY) 3 4 #remove small dots 33
Figure 29: Histogram for Figure 28. 5 thresh1 = cv2.dilate(thresh1, np.ones((5, 5), np.uint8)) 6 thresh1 = cv2.erode(thresh1, np.ones((5, 5), np.uint8)) 7 8 cv2.imwrite("output/Example1BWThresholding.png", thresh1) that produce Figure 30. We are ready to detect contours and draw minimum area rectangle 1 #reverse colors in mask 2 thresh2 = abs(thresh1 - 255) 3 4 #find contours 5 contours, _ = cv2.findContours(thresh2, cv2.RETR_LIST, 6 cv2.CHAIN_APPROX_TC89_L1 ) 7 blank = np.ones(thresh2.shape)*255 8 cv2.drawContours(blank, contours, -1, (0, 0, 0), 3) 9 print("# of contours found: ", len(contours)) 10 cv2.imwrite("output/Example1BWCountour.png", blank) 11 12 #minmimal area rectangle 13 img1_tmp1 = img1.copy() 14 img1_tmp2 = img1.copy() 15 for idx, cnt in enumerate(contours): 16 minRect = cv2.minAreaRect(cnt) 17 box = cv2.boxPoints(minRect) 34
Figure 30: The result of thresholding in Figure 28. 18 box = np.int0(box) 19 cv2.drawContours(img1_tmp1, [box], 0, (0,0,255), 2) 20 cv2.drawContours(img1_tmp2, [box], 0, (0,0,255), 2) 21 cv2.putText(img1_tmp2, '{}'.format(idx), (int(box[0][0]), 22 int(box[0][1]) ), cv2.FONT_HERSHEY_SIMPLEX, 0.5, 23 (255, 255, 255), 1, cv2.LINE_AA) 24 cv2.imwrite("output/Example1BWMinAreaRectangle.png", img1_tmp1) 25 cv2.imwrite("output/Example1BWMinAreaRectangleWithIndices.png", 26 img1_tmp2) That produces Figure 31. We can see that most of the granules were detected by this (semi) automatic procedure. The final step is to fit ellipses to the contours, and save their parameters in CSV file, that is realized by the following code 1 blank = np.ones(thresh2.shape)*255 2 ellipsesData = [] 3 img_tmp = img1.copy() 4 for cnt in contours: 5 #if less points preventing fitting ellipse 6 if cnt.size < 10 : 7 continue 8 x,y,w,h = cv2.boundingRect(cnt) 9 #remove small conoturs 10 if w < 2 or h < 2: 35
Figure 31: Minimal area rectangles enclosing detected contours with enu- meration. 11 continue 12 ellipse = cv2.fitEllipse(cnt) 13 cv2.ellipse(blank, ellipse, (0,255, 255), 1, cv2.LINE_AA) 14 cv2.ellipse(img_tmp, ellipse, (0,255, 255), 1, cv2.LINE_AA) 15 ellipsesData.append({"id": idx, "x":ellipse[0][0], 16 "y":ellipse[0][1], "majorAxis": ellipse[1][0], "minorAxis": 17 ellipse[1][1], "angle": ellipse[2]}) 18 19 cv2.imwrite("output/Example1BWEllipses.png", blank) 20 cv2.imwrite("output/Example1BWEllipsesColor.png", img_tmp) 21 22 df_ellipses = pd.DataFrame.from_records(ellipsesData) 23 df_ellipses.to_csv("output/ellipses.csv") 24 25 plt.hist(np.array(df_ellipses.angle)%90) 26 plt.xlabel("angle[deg] mod 90") 27 plt.ylabel("No. of ellipses") 28 plt.savefig("output/Example1AngleDistibution.png") 29 plt.show() This produces Figure 32. Moreover, when we are looking for the orientation of ellipses that can represent the orientation of non-spherical granules, we get from the above code Figure 33. 36
Figure 32: Ellipses fit to the contours. Figure 33: Angle distributions for ellipses. 37
2.2 Example 2 - classification using images We now want to classify contours from the previous example using the k- means algorithm. As mentioned above, the classification results must be consulted with an expert. The first part of the code related to preparing images and creating masks is the same as before. We provide it for com- pleteness 1 import numpy as np 2 import pandas as pd 3 import matplotlib.pylab as plt 4 import cv2 5 6 img1 = cv2.imread("Picture1.png") 7 8 #cut out printout 9 img1 = img1[50:-10,10:] 10 cv2.imwrite("output/Example2Cut.png", img1) 11 12 #histograms of colors 13 hist= cv2.calcHist([img1],[0],None,[256],[0,256]) 14 hist2=cv2.calcHist([img1],[1],None,[256],[0,256]) 15 hist3=cv2.calcHist([img1],[2],None,[256],[0,256]) 16 ##plot histograms 17 plt.subplot(1, 3, 1) 18 plt.plot(hist/max(hist)) 19 plt.title('Hist Blue : ') 20 plt.subplot(1, 3, 2) 21 plt.plot(hist2/max(hist2),color="green") 22 plt.title('Hist Green : ') 23 plt.subplot(1, 3, 3) 24 plt.plot(hist3/max(hist3),color="red") 25 plt.title('Hist Red : ') 26 plt.savefig("output/Example2ColorHistograms") 27 plt.show() 28 29 #histogram equalization 30 src = cv2.cvtColor(img1, cv2.COLOR_BGR2HSV) 31 h,s,v = cv2.split(src) 32 v = cv2.equalizeHist(v) 33 src = cv2.merge([h,s,v]) 34 img1 = cv2.cvtColor(src, cv2.COLOR_HSV2BGR) 38
35 cv2.imwrite("output/Example2Equalization.png", img1) 36 37 #histograms of colors after normalization 38 hist= cv2.calcHist([img1],[0],None,[256],[0,256]) 39 hist2=cv2.calcHist([img1],[1],None,[256],[0,256]) 40 hist3=cv2.calcHist([img1],[2],None,[256],[0,256]) 41 plt.subplot(1, 3, 1) 42 plt.plot(hist/max(hist)) 43 plt.title('Hist Blue : ') 44 plt.subplot(1, 3, 2) 45 plt.plot(hist2/max(hist2),color="green") 46 plt.title('Hist Green : ') 47 plt.subplot(1, 3, 3) 48 plt.plot(hist3/max(hist3),color="red") 49 plt.title('Hist Red : ') 50 plt.savefig("output/Example2ColorHistogramsEqualization.png") 51 plt.show() 52 53 #scale figure 54 img1 = cv2.resize(img1,None,fx=4, fy=4, interpolation = 55 cv2.INTER_CUBIC) 56 # convert to BW 57 img1_BW = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY) 58 cv2.imwrite("output/Example2BW.png", img1_BW) 59 60 #bilateral filtering 61 img1_BW = cv2.bilateralFilter(img1_BW,9,75,75) 62 cv2.imwrite("output/Example2BWFiltering.png", img1_BW) 63 64 #compute gray histogram 65 hist1 = cv2.calcHist([img1_BW], [0], None, [256], [0,256]) 66 plt.plot(hist1.ravel()/hist1.max()) 67 plt.savefig("output/Example2HistogramBW.png") 68 plt.show() 69 70 #tresholding - we select 40 for selecting black granules 71 ret1m, thresh1 = cv2.threshold(img1_BW, 40, 256, 72 cv2.THRESH_BINARY) 73 74 #remove small dots 75 thresh1 = cv2.dilate(thresh1, np.ones((5, 5), np.uint8)) 39
76 thresh1 = cv2.erode(thresh1, np.ones((5, 5), np.uint8)) 77 78 cv2.imwrite("output/Example2BWThresholding.png", thresh1) 79 80 #reverse colors in mask 81 thresh2 = abs(thresh1 - 255) 82 83 #find contours 84 contours, _ = cv2.findContours(thresh2, cv2.RETR_LIST, 85 cv2.CHAIN_APPROX_TC89_L1 ) 86 blank = np.ones(thresh2.shape)*255 87 cv2.drawContours(blank, contours, -1, (0, 0, 0), 3) 88 print("# of contours found: ", len(contours)) 89 cv2.imwrite("output/Example2BWCountour.png", blank) 90 91 #minmimal area rectangle 92 img1_tmp1 = img1.copy() 93 img1_tmp2 = img1.copy() 94 for idx, cnt in enumerate(contours): 95 minRect = cv2.minAreaRect(cnt) 96 box = cv2.boxPoints(minRect) 97 box = np.int0(box) 98 cv2.drawContours(img1_tmp1, [box], 0, (0,0,255), 2) 99 cv2.drawContours(img1_tmp2, [box], 0, (0,0,255), 2) 100 cv2.putText(img1_tmp2, '{}'.format(idx), (int(box[0][0]), 101 int(box[0][1]) ), cv2.FONT_HERSHEY_SIMPLEX, 0.5, 102 (255, 255, 255), 1, cv2.LINE_AA) 103 cv2.imwrite("output/Example2BWMinAreaRectangle.png", img1_tmp1) 104 cv2.imwrite("output/Example2BWMinAreaRectangleWithIndices.png", 105 img1_tmp2) We now crop from the initial image as bounding rectangles surrounding contours. These images will be reshaped into vectors to be an input for the ML algorithm. The following code prepares the data 1 image = img1.copy() 2 img_tmp = img1.copy() 3 filtered_contours = [] 4 selected_img = [] 5 for idx, contour in enumerate(contours): 6 #draw selected contour 40
7 minRect = cv2.minAreaRect(contour) 8 box = cv2.boxPoints(minRect) 9 box = np.int0(box) 10 cv2.drawContours(img_tmp, [box], -1, (0,0,255), 1) 11 12 #extract image from bounding rectangle 13 x, y, w, h = cv2.boundingRect(contour) 14 cropped = image[y:y + h, x:x + w, :].copy() 15 feature = cv2.resize(cropped, (20, 20)) 16 cv2.imwrite("output/figs/{}.png".format(idx), feature) 17 18 #scale cropped image to 1D array 19 feature = feature.reshape(-1) 20 selected_img.append(feature) 21 cv2.imwrite("output/Exmaple2SelectedContours.png", img_tmp) 22 df = pd.DataFrame(selected_img) We now use these data to determine optimal number of clusters (similar granules) 1 from sklearn.cluster import KMeans 2 k = range(1,20) 3 WSSE = [] 4 for i in k: 5 model = KMeans(n_clusters = i, n_init='auto') 6 model.fit_predict(df) 7 WSSE.append(model.inertia_) 8 plt.plot(k,WSSE) 9 plt.xlabel('K') 10 plt.ylabel('WSSE') 11 plt.grid(True) 12 plt.savefig("output/Example2kmeans.png") 13 plt.show() that produces Figure 34. From the ’elbow curve’, we can observe that k = 4 clusters are sufficient. We therefore do classification using this value. More- over, we fit ellipses to the contours and assign them to the class labels. The following lengthy code does the mentioned operations 1 km = KMeans( n_clusters=4, n_init='auto') 2 df['label'] = km.fit_predict(df) 41
Figure 34: Angle distributions for ellipses. 3 4 df.to_csv("output/LabelledContours.csv") 5 6 img_tmp = image.copy() 7 img_ellipse = image.copy() 8 img_ellipse_closed = image.copy() 9 ellipsesData = [] 10 for label, df_grouped in df.groupby('label'): 11 Ngranules = len(df_grouped) 12 print(Ngranules, " in label ", label) 13 mask = np.zeros_like(image[:, :, 0]) 14 cv2.drawContours(mask, [contours[i] for i in 15 df_grouped.index], -1, (255), -1) 16 masked_image = cv2.bitwise_and(image, image, mask=mask) 17 masked_image[mask==0] = (255,255,255) 18 cv2.imwrite("output/Example2label{}.png".format(label), 19 masked_image) 20 21 print("fitting for label = ", label) 22 for idx in df_grouped.index: 23 minRect = cv2.minAreaRect(contours[idx]) 24 box = cv2.boxPoints(minRect) 25 box = np.int0(box) 42
26 cv2.drawContours(img_tmp, contours[idx], -1, 27 (0,0,label*50), 1) 28 29 #fit ellypses 30 if contours[idx].size < 10 : 31 continue 32 if w < 2 or h < 2: 33 continue 34 ellipse = cv2.fitEllipse(contours[idx]) 35 ellipsesData.append({"id": idx, "x":ellipse[0][0], 36 "y":ellipse[0][1], "majorAxis": ellipse[1][0], 37 "minorAxis": ellipse[1][1], "angle": ellipse[2], 38 "class": label}) 39 if max(ellipse[1]) < 50: 40 cv2.ellipse(img_ellipse, ellipse, (0,0,label*50), 41 1, cv2.LINE_AA) 42 cv2.putText(img_ellipse, '{}_{}'.format(idx, label), 43 (int(box[0][0]), int(box[0][1]) ), 44 cv2.FONT_HERSHEY_SIMPLEX, 0.5, (255, 255, 255), 1, 45 cv2.LINE_AA) 46 cv2.ellipse(img_ellipse_closed, ellipse, 47 (0,0,label*50), -1, cv2.LINE_AA) #open ellipses 48 cv2.putText(img_ellipse_closed, '{}_{}'.format(idx, 49 label), (int(box[0][0]), int(box[0][1]) ), 50 cv2.FONT_HERSHEY_SIMPLEX, 0.5, (255, 255, 255), 1, 51 cv2.LINE_AA) 52 53 df_ellipses = pd.DataFrame.from_records(ellipsesData) 54 df_ellipses.to_csv("output/ellipses.csv") 55 56 cv2.imwrite('output/Example2_color_segmentation.png', img_tmp) 57 cv2.imwrite('output/Example2_color_ellypses.png', img_ellipse) 58 cv2.imwrite('output/Example2_color_ellypses_closed.png', 59 img_ellipse_closed) which produces figures marking classified ellipses, as Figure 35. 2.3 Example 3 In this example, we do a similar analysis to the previous one using Figure 36. After reading and cropping the image by the following code snippet 43
Figure 35: Fitted ellipses with numeration and class in the format: num- ber label. Figure 36: Initial image. Microscope image of cancerous mouse brain biopsy. 44
Figure 37: Cropped image. 1 img1 = img1[50:-80,10:] that result in Figure 37. The histogram of intensity for a grayscale image, Figure 38 suggests the threshold at about 50. The identified contours are presented in Figure 39. Finally, on classifying and fitting ellipse we get Figure 40. 2.4 Example 4 - classification using mean color and moments of contours We use again the image from Figure 36 and the code from the previous ex- ample. However, instead of using a full image representing a granule cropped from the original image as data for the classifier, we use mean color (each color channel) and Hu moment invariants. These parameters can be impor- tant if we expect that the geometry of granules, as well as their (average) color, are parameters that are important to discriminating different classes. The change in the previous code that prepares the data is as follows 1 import math 2 image = img1.copy() 3 img_tmp = img1.copy() 4 filtered_contours = [] 45
Figure 38: Intensity histogram. Figure 39: Contours marked by minimal area rectangles and enumerated. 46
Figure 40: Ellipses fitted to contours with classes marked by different colors and enumeration in the format: number label. 5 selected_data = [] 6 for idx, contour in enumerate(contours): 7 # get mean color of contour: 8 masked = np.zeros_like(image[:, :, 0]) 9 cv2.drawContours(masked, [contour], 0, 255, -1) 10 11 B_mean, G_mean, R_mean, _ = cv2.mean(image, mask=masked) 12 13 #calculate moments 14 moments = cv2.moments(contour) 15 # Calculate Hu Moments 16 huMoments = cv2.HuMoments(moments) 17 18 #scaling Hu moments 19 for i in range(0,7): 20 #print("hu moment ", i, " ", huMoments[i]) 21 if huMoments[i] == 0: 22 huMoments[i] = 0.0 23 else: 24 huMoments[i] = -1* math.copysign(1.0, huMoments[i]) 25 * math.log10(abs(huMoments[i])) 26 27 selected_data.append({'B_mean': B_mean, 'G_mean': G_mean, 47
Figure 41: Ellipses fitted to contours with classes marked by different colors and enumeration in the format: number label. 28 'R_mean': R_mean, 'h0':huMoments[0], 'h1':huMoments[1], 29 'h2':huMoments[2], 'h3':huMoments[3], 'h4':huMoments[4], 30 'h5':huMoments[5], 'h6':huMoments[6]}) 31 cv2.imwrite("output/Exmaple2SelectedContours.png", img_tmp) 32 df = pd.DataFrame(selected_data) For comparison, fitted ellipses are presented in Figure 41. 3 Automatic detection for zones in tissue The structural sums serve as the cornerstone for mathematical models of random structures within the framework of the analytical Representative Volume Element (aRVE) theory. We quantified the inter-phase interactions within media and derived analytical formulas for their macroscopic properties based on the Riemann-Hilbert and R-linear problems for analytic functions. 48
Figure 42: Cancerous tissues. The picture is due to Pablo Losa represented in the COST Action CA23125 workshop in Brcelona 2026. 49
Figure 43: The processing of the central square in Figure 42. White areas and their gravitational centers (red dots). The main problem of homogenization (global/local) is to replace a hetero- geneous medium with an equivalent homogeneous medium and to determine its overall behavior. Following Parts 2.1 and 2.2 consider a doubly periodic multi-phase com- posite when the host G+m1ω1 +m2ω2 and the inclusions Gk +m1ω1 +m2ω2 are occupied by dielectric materials. Let the permittivity of the host be nor- malized to unity, and the permittivity of the kth inclusion be a complex number εk. Let u and uk denote the potentials in G and Gk, respectively. The complex functions u and uk satisfy Laplace’s equation in the correspond- ing domains and are continuously differentiable in their closures. The perfect contact between the components is expressed by the boundary condition u(t) = uk(t), ∂u ∂n(t) = εk ∂uk ∂n (t), t ∈ ∂Gk (k = 1, 2, . . . , N ). (3.1) 50
Following Parts 2.1 and 2.2 consider the contrast parameter ϱk = εk − 1 εk + 1 . (3.2) Introduce the averaged contrast parameter over inclusions ⟨ϱ⟩ = NX k=1 ϱk|Gk|. Let ak be the complex coordinate of the gravity center of Gk; E2(z) = ℘(z) + π the Eisenstein function; Emk := E2(am − ak) for m̸ = k. Introduce L = 1 π NX k,m=1 ϱkϱm|Gk||Gm|Emk. (3.3) Introduce the shape cumulative parameter J = NX k=1 ϱ2 kJk, Jk = 1 4π Z Lk dt Z Lk τ τ − t dτ. (3.4) The above sums can be written in dimensionless form using the total fraction f of inclusions J0 = f −1J , L0 = f −2L. The permittivity tensor can be calculated by the formula ε⊥ = (1 + 2⟨ϱ⟩)I + 2 PN k=1 ϱ2 k Re Jk −Im Jk −Im Jk −Re Jk + 2 π PN k,m=1 ϱkϱm|Gk||Gm| Re Emk − Im Emk −Im Emk 2 − Re Emk + O(ϱ3 0f 3), (3.5) where ϱ0 = maxk |ϱk|. We us the structural sums em1,...,mq = 1 N 1+ 1 2 (m1+···+mq ) × P k0,...,kq Em1 (ak0 − ak1 )CEm2 (ak1 − ak2 ) . . . Cq+1Emq (akq−1 − akq ), (3.6) where Ep(z) are the Eisenstein functions, C the complex conjugation. For instance, e2 and e22 take the following form e2 = 1 N 2 NX k0=1 NX k1=1 E2(ak0 − ak1 ), e22 = 1 N 3 NX k0=1 NX k1=1 NX k2=1 E2(ak0 − ak1 )E2(ak1 − ak2 ). (3.7) 51
The structural sum (3.6) is a discrete multiple convolution of the Eisenstein functions. A fast, almost linear algorithm in N for its computation was developed in [6, 7]. The theory of aRVE can be considered as a constructive application of the Decomposition Theorem. One of aRVE applications is outlined below. Let us have at our disposal two plane pictures of dispersed composites and their digital treatment in the form of two sets of the centers of inclusions a = {a1, a2, . . . , aN } and a′ = {a′ 1, a′ 2, . . . , a′ N ′ } scaled to the periodic unit cell. We want to know whether these two composites belong to the same class of random composites. First, we must calculate all the structural sums em1,...,ms for two sets of points a and a′. The coefficient Ak on f k+1 is a linear combination of the structural sums em1,...,ms with m1 + . . . + ms = 2k. Therefore, we have to compare the set of complex values E = {E1, E2, . . .} where E1 = {e2}, E2 = {e22}, E3 = {e33, e222}, E4 = {e44, e332, e233, e2222}, . . . . (3.8) constructed for a and the set E′ for a′. The set E is infinite and may be truncated by the concentration tolerance. Let it be equal to O(f 8) for definiteness. Then, we have to calculate the finite set of structural sums em1,...,ms E1 = {e2}, E2 = {e22}, . . . , E7 = {e77, e266, . . . , e2222222}. (3.9) It is worth noting that the 2-point correlation functions yield the same result as only one constant e2; the 3-point correlation functions yield the the same result as the infinitely many elements epp (p = 2, 3, . . .); the 4- point correlation functions yield the same result as the infinitely many triple elements em1,m2,m3 from E and so forth. Using aRVE is equivalent to the usage of the multi-point correlation functions. The structural sums for the random i.i.d. simulations e′ 2 = 3.14, e′ 22 = 11.5249, e′ 33 = −23.3673, e′ 44 = 7.9885. (3.10) The structural sums for the considered structure has the form e2 = 3.3876 − 0.3697i, e22 = 41.335, e33 = −338.484, e44 = 6541.22. (3.11) 4 Conclusion and practical application The theory of an aRVE bridges the gap between classical homogenization theory and the incorporation of randomness in heterogeneous media. Its 52
methodological and practical application enables a systematic treatment of the variability in random composites within the framework of homogenization principles and asymptotic analysis. Application to the diagnosis and classification of glioma structures re- quires a set of corresponding pictures to apply the developed machine learn- ing classification algorithm. The following strategy of cancer detection can be applied by by the de- veloped algorithms to obtain practical results: • images of cancer samples with the precise description of zones with their permittivity data; • processing the pictures by the presented algorithm; • application of math formulas to permittivity and structural sums for a classification. In order to perform the above strategy step we need: • physical constants (perhaps as functions of frequency) of the image components (2 or multi-phase). It must be precisely described (by color) or by direct indication, which domains have the given constant. • to present at least 10 pictures of healthy and 10 pictures of diseased samples. 10 + 10 pictures should present one type of healthy - dis- eased samples to process their and use in the next classification of new pictures by supervised machine learning. References [1] Example code from this book:. https://github.com/rkycia/ AnisotropyOfMetamaterialsBook. Accessed: 27.09.2024. [2] Opencv library. https://opencv.org/. Accessed: 29.07.2024. [3] Scientific python lectures. https://lectures.scientific-python. org/. Accessed: 29.07.2024. [4] Scikit learn library. https://scikit-learn.org. Accessed: 29.07.2024. [5] Uci ml hand-written digits datasets. https://archive.ics.uci.edu/ dataset/80/optical+recognition+of+handwritten+digits. Ac- cessed: 29.07.2024. 53
[6] R. Czapla, W. Nawalaniec, and V. Mityushev. Effective conductivity of random two-dimensional composites with circular non-overlapping inclusions. Computational Materials Science, 63:118–126, 2012. [7] P. Dryga´s, S. Gluzman, V. Mityushev, and W. Nawalaniec. Applied analysis of composite media: analytical and computational results for materials scientists and engineers. Woodhead Publishing, 2019. [8] R.O. Duda and P.E. Hart. Use of the hough transformation to detect lines and curves in pictures. Communications of the ACM, 15(1):11–15, 1972. [9] J. Flusser. On the independence of rotation moment invariants. Pattern recognition, 33(9):1405–1410, 2000. [10] J. Flusser and T. Suk. Rotation moment invariants for recogni- tion of symmetric objects. IEEE Transactions on Image Processing, 15(12):3784–3790, 2006. [11] Aur´elien G´eron. Hands-on machine learning with Scikit-Learn, Keras, and TensorFlow. ” O’Reilly Media, Inc.”, 2022. [12] R.C. Gonzalez. Digital image processing. Pearson education india, 2009. [13] Joseph Howse and Joe Minichino. Learning OpenCV 4 Computer Vision with Python 3: Get to grips with tools, techniques, and algorithms for computer vision and machine learning. Packt Publishing Ltd, 2020. [14] Ming-Kuei Hu. Visual pattern recognition by moment invariants. IRE transactions on information theory, 8(2):179–187, 1962. [15] A. Kaehler and G. Bradski. Learning OpenCV 3: computer vision in C++ with the OpenCV library. ” O’Reilly Media, Inc.”, 2016. [16] R. Klette. Concise computer vision, volume 233. Springer, 2014. [17] V. Mityushev, Kycia R., W. Nawalaniec, and N. Rylko. Introduction to Mathematical Modeling and Computer Simulations. Francis & Taylor, 2nd edition, 2025. [18] N. Otsu et al. A threshold selection method from gray-level histograms. Automatica, 11(285-296):23–27, 1975. [19] S. Raschka, Yuxi Hayden Liu, and V. Mirjalili. Machine Learning with PyTorch and Scikit-Learn: Develop machine learning and deep learning models with Python. Packt Publishing Ltd, 2022. 54
[20] C.-H. Teh and R.T. Chin. On the detection of dominant points on digital curves. IEEE Transactions on pattern analysis and machine intelligence, 11(8):859–872, 1989. [21] A.F. Vill´an. Mastering OpenCV 4 with Python: a practical guide cov- ering topics from image processing, augmented reality to deep learning with OpenCV 4 and Python 3.7. Packt Publishing Ltd, 2019. 55