Development and testing of a novel image analysis algorithm for descriptive evaluation of shape change of a shrinkable soft material

Soft material can undergo non-uniform deformation or change of shape upon processing. Identifying shape and its change is nevertheless not straightforward. In this study, novel image-based algorithm that can be used to identify shapes of input images and at the same time classify non-uniform deformation into various patterns, i.e., swelling/shrinkage, horizontal and vertical elongations/contractions as well as convexity and concavity, is proposed. The algorithm was first tested with computer-generated images and later applied to agar cubes, which were used as model shrinkable soft material, undergoing drying at different temperatures. Shape parameters and shape-parameter based algorithm as well as convolutional neural networks (CNNs) either incorrectly identified some complicated shapes or could only identify the point where non-uniform deformation started to take place; CNNs lacked ability to describe non-uniform deformation evolution. Shape identification accuracy of the newly developed algorithm against computer-generated images was 65.88%, while those of the other tested algorithms ranged from 34.76 to 97.88%. However, when being applied to the deformation of agar cubes, the developed algorithm performed superiorly to the others. The proposed algorithm could both identify the shapes and describe their changes. The interpretation agreed well with that via visual observation.


Results and discussion
Test of shape identification algorithms. First, the accuracies of the newly proposed shape identification procedures with and without image alignment were compared; the results are shown in Table 1. The results of the use of the algorithm by Igathinathane et al. 17 and CNNs are also listed. Note that the shape identification accuracy was calculated by comparing the known shapes of the computer-generated images with those identified by each tested algorithm. A count was made on both correct and incorrect shape identifications; percent accuracy was then calculated. In the case of the newly developed algorithm, identifying the shape of input images without any alignment yielded the lowest accuracy, as expected. Unaligned images had higher chances of being overlapped with other reference geometries than with the most suitable ones. Nevertheless, although the alignment step could help improve the accuracy of shape identification, such a step resulted in a longer processing time. Zigzagging artefact, which was the results of the alignment step 28 , might also affect shape identification. Different alignment methods (see Fig. 1) yielded different levels of accuracy as shown in Table 1. Equivalent ellipsis method resulted in the lowest accuracy. When images were not deformed (0, 0, 0, 0, 0), major axis of those images could have already been aligned, either on the horizontal or vertical axis, causing the method not to rotate the images until reaching α = 0° or even rotating the images upside-down; this is because α is always the acute angle between the major axis and horizontal axis. When applying the equivalent ellipsis method to images with horizontal elongation (+, +, 0, 0, 1), the method yielded the lowest accuracy. For example, when a triangle or pentagon was horizontally elongated, the method identified it instead as a rectangle or hexagon. Such an error occurred because the apex angle would change from smaller than 90° to larger than 90°. The method performed better when images were vertically elongated (+, 0, +, 0, 1) because the major axis of those images would align on the longest axis of the images. Nevertheless, this alignment method required the shortest processing time (around 0.25 s per image using the currently employed hardware) because it refers to only one angle (α) in order to perform the alignment process.
Extrema method was noted to give more accurate identification results than the equivalent ellipsis method. The accuracy was higher than 50% if input images resembled their reference geometric shapes. This is because the extrema points of the images would be the same as those of the reference geometric images. On the other hand, when images were convexed (0, 0, 0, +, 1), the extrema points sometimes changed to the convex points. Convexity also caused images to look more like a circle or an ellipse. This alignment method required longer processing time than the equivalent ellipsis method (around threefold longer) as the former needs to refer to up to four angles to perform the alignment process (see Fig. 2).
Manual alignment method gave the highest accuracy (higher than 60% in all cases). The major advantage of this method arose when images experienced horizontal (+, +, 0, 0, 1) or vertical elongation (+, 0, +, 0, 1); such images would suffer less inaccuracy than the ones with no deformation (0, 0, 0, 0, 0). This is because the manual alignment method rotates an image from 0° to 360° at 1° interval, causing some zigzagging artefact on the image edge. However, when images were elongated, the dimension ratio of the artefact to the characteristic length of the images would decrease. Nevertheless, this method required the longest time, 51.29 s per image or around 70-fold longer than the extrema method. This is expected as it has 360 possible angles to refer to.
Although the manual alignment method resulted in the highest accuracy, the processing time was excessive. This is clearly not appropriate, especially if real-time image-based process control is to be conducted. For this reason, extrema method was chosen because it exhibited higher accuracy than the equivalent ellipsis method and yet performed 70 times faster than the manual alignment method.
It is seen in Table 1 that the developed algorithm (employing the extrema alignment method) gave 100.00% and 92.92% accuracy when input images did not suffer any deformation (0, 0, 0, 0, 0) and when images were concaved (0, 0, 0, −, 1), respectively. This is because the points at the angles of those images were distinct (except in the case of circles). This in turn caused the RMSD values to be lower. Swollen images (+, , +, 0, 0) having more zigzagging edges, caused by the forced deformation process, would, for example, result in pentagons or hexagons  www.nature.com/scientificreports/ looking more like circles. Horizontal elongation exhibited stronger adverse effect than vertical elongation because when images with apex angles of less than 90° (such as triangles and pentagons) were horizontally elongated, the angles would become larger than 90°. This in turn affected RMSD because this value was calculated from d values at the same angle θ. The algorithm did not well identify images when they were convexed as the image would become more like circles. When all 4200 computer-generated images were analyzed using the developed algorithm with the extrema method, the overall capability to accurately identify the shapes of input images was noted to be 65.88%. In contrast, when the algorithm by Igathinathane et al. 17 and CNN with ResNet-50 architecture (CNN ResNet-50) and CNN with Xception architecture (CNN Xception) were used for the same shape identification purpose, the overall accuracies were noted to be 34.76%, 85.48 and 97.88% respectively. The reason for the much lower overall accuracy of the algorithm by Igathinathane et al. 17 is that the algorithm could not identify pentagons and hexagons as these shapes are not included in the algorithm; computer-generated pentagons and hexagons were therefore incorrectly identified as either ellipses or rectangles. This clearly indicates the disadvantage of this (or similar algorithms) where new reference geometric shapes always need to be added and analyzed when dealing with input images of peculiar shapes. In the case of CNNs, the algorithms were not good at identifying both horizontally and vertically convexed images. This may probably be due to the inherently ambiguous sorting classes and blurring appearance 29 . CNN Xception performed better than CNN ResNet-50, as expected 30 . It is important to note, however, that both the newly developed algorithm and that of Igathinathane et al. 17 require no training, while CNN ResNet-50 and CNN Xception require significant training time. Using the currently employed hardware, as much as 41,942 and 138,269 s were required, respectively, for the two CNNs.
Comparative performance of developed algorithm and existing shape parameters. Effect of size change. The effect of the change in image size on the Extent was first evaluated. Size change did not significantly affect the Extent. This is because when the size of an image increased, the size of the rectangular bounding box also increased at the same rate. Each geometric shape was noted to exhibit its own Extent value; the values for a triangle, rectangle, pentagon, hexagon and circle are 0.525 ± 0.003, 1.000 ± 0.000, 0.710 ± 0.002, 0.768 ± 0.002 and 0.805 ± 0.004, respectively. Extent nevertheless has a significant disadvantage that it can only be used to identify reference geometric shapes. This shape parameter cannot be used to describe non-uniform deformation; only a point where non-uniform deformation starts to take place could be identified 15 .
In the case of fractal dimension, change of this parameter could be divided into three periods based on the ratio of area of interest (AOI) to the total area of an image. When the aforementioned ratio was lower than 5%, fractal dimension was close to unity. This unexpected result was noted because the method that was used to calculate the fractal dimension here was the box counting method. When the size decreased, AOI resembled a point. This implies that fractal dimension calculated by the box counting method could not be used when the ratio of AOI to the total area of an image is lower than 5%. On the other hand, when such a ratio was higher than 5% but lower than 20%, fractal dimensions of a triangle, rectangle, pentagon, hexagon and circle were noted to be 1.826 ± 0.061, 1.863 ± 0.036, 1.852 ± 0.052, 1.869 ± 0.051, 1.883 ± 0.051, respectively. Finally, when such a ratio was higher than 20%, all fractal dimensions approached 2.
In the case of the presently developed algorithm, the change in size affected i, j and k. This is expected as i is related to image area, which is in turn related to size; j and k are also linearly related to horizontal and vertical dimension changes. Parameter l, which is used to indicate convexity and concavity was not affected by the change in size, as expected.
Effects of horizontal and vertical length changes. Both horizontal elongation and contraction did not affect the Extent because these changes equally affected both AOI and its rectangular bounding box. In the case of the newly developed algorithm, horizontal elongation and contraction affected i and j but did not affect k and l. When an image was horizontally elongated or contracted, its size also increased or decreased, respectively. However, since the change occurred only in the horizontal direction, k was not affected. Parameter l was also not affected as no convexity took place. The above-mentioned reasons can also be used to explain the effect of vertical length change.
Effects of convexity and concavity. Extent of a triangle was most significantly affected and so can be used to quantify both the convexity and concavity. In contrast, the Extent of a rectangle exhibited similar patterns of changes, whether convexity or concavity took place; Extent cannot therefore be used to classify convexity and concavity of this geometric shape. Extent values of a concaved pentagon, hexagon and circle exhibited similar patterns; however, the values belonging to these shapes that were made to convex remained almost unchanged. Therefore, the Extent can only be used to identify concavity of these shapes.
Since each geometric shape exhibited its own Extent value, convexity and concavity of a triangle could, for example, be identified by looking at the deviation of its Extent value from the reference value of 0.525; higher Extent value than the reference value indicated convexity and vice versa. Similar identifications can also be made for a pentagon, hexagon and circle. Since convexity and concavity caused the image perimeter to be rougher, fractal dimension increased. However, this parameter cannot be used to describe such a change of shape, only being able to identify if any change exists.
In the case of the presently developed algorithm, concavity and convexity only affected parameter l. As a result, this parameter can be used to identify and quantify convexity and concavity by looking at the deviation of its value from zero in much the same way as in the case of Extent. However, this parameter can be applied to all tested shapes, including a rectangle. www.nature.com/scientificreports/ Note that parameter m is used to identify if deformation is uniform or non-uniform. Such a decision is made based on the values of α and l. Uniform deformation is said to take place only if both the value of α is higher than 95% and the value of l is higher than 99%. The chosen threshold of 95% is based on the typical confidence level used in most engineering statistical analysis. The threshold of 99% was, on the other hand, chosen by adjusting the value until the computed identification matched that visually observed.
Test of developed algorithm on agar cubes. Performance of the developed algorithm was evaluated by applying it to agar cubes containing either 0 or 20% sucrose undergoing hot air drying at 80 °C (see Fig. 3). Top-view images of agar cubes containing 0 and 20% sucrose are shown in Table 2. The algorithm indicated that the top-view shape of agar cubes with 0% sucrose changed from square to rectangle to hexagon when the moisture ratios were about 0.78 and 0.36, while no shape change was detected in the case of agar cubes with 20% sucrose. Agar with 0% sucrose had almost no solids to support its structure after water had been evaporated. On the other hand, agar with 20% sucrose possessed more solids to support its structure and hence deformed more uniformly. In  www.nature.com/scientificreports/ other words, the cubes with 20% sucrose experienced only the change in size; no convexity and concavity were observed. Top-view shape changes of agar cubes with 0% sucrose could be divided into three periods. First, the shape of the cubes was slightly convexed. This is because their edges were hardened while the other portions, which were still moist and elastic, were forced to convexed out (see the results at the moisture ratios of 0.78 and 0.56). i decreased from 100 to 84.35 and 68.69% when the area decreased by 15.65 and 31.31%, while j decreased from 100 to 96.84 and 84.98%, indicating the horizontal contractions of 3.16 and 15.02%. k also decreased from 100 to 96.93 and 81.91%, indicating the vertical contractions of 3.07 and 18.09%. The decreasing rate of k was higher than that of j, implying that the vertical contraction was more extensive than the horizontal one. This is because the direction of drying air flow in the present experimental study was perpendicular to the top horizontal edges of the cubes. In the case of l, its value increased from 100 to 101.11 and 102.76%, indicating, as mentioned earlier, that convexity had occurred. The top-view shape of the cubes then appeared as rectangle rather than as square. Parameter m also correctly indicated the non-uniform deformation.
During the second period, when the moisture ratios were between 0.56 and 0.36, the shape of the cubes changed into hexagon. i decreased from 68.69 to 52.75%, indicating that shrinkage still occurred. The value of j decreased from 84.98 to 77.47%, while that of k decreased from 83.91 to 75.48%. Trends of changes of j and k were similar to those during the first period in that there was more extensive vertical contraction than the horizontal one. The value of l changed from 102.76 to 103.32%, indicating that there was also convexity. It was indeed this convexity that led to further shape change. Parameter m again indicated non-uniform deformation.
Finally, when the moisture ratios were lower than 0.36, structural collapse was noted, resulting in the change of shape of agar cubes from hexagon to pentagon. This is because agar cubes exhibited concavity in the middle softest region prior to finally being hardened (see the results at the moisture ratios of 0.20 and 0.06). The value of i decreased to 24.07%, implying that agar cubes extensively shrunk; j and k also decreased to 58.01 and 50.19%, respectively. The value of l decreased to 99.33%, implying that the cubes were concaved. Parameter m again indicated that the deformation was non-uniform in nature.
The results obtained by the developed algorithm were compared with those obtained using the algorithm developed by Igathinathane et al. 17 and CNN Xception, in terms of Shape Ig and Shape CNN , respectively; CNN Xception was selected due to its higher overall identification accuracy. The results obtained using Extent and fractal dimension were also compared. The results are as shown in Table 2. The developed algorithm could both identify and describe the changes of shape. The interpretation indeed agreed with that via visual observation. On the contrary, the algorithm of Igathinathane et al. 17 could only identify the shape as "inclined rectangle" even though the cubes changed into other shapes (e.g., pentagon, hexagon) rather than rectangle. Interestingly, although CNN Xception could identify the shapes of computer-generated images very accurately, it could only identify the shape of agar cubes as "rectangle. " This might be due to the insufficient number of images for the network training. Although as many as 45,000 images were used for training the network, it still did not succeed in this real-world application. This illustrates a limitation of CNNs when being applied even to identify simple shapes of practical materials and/or objects. Acquiring a sufficient number of real-world images for training may not always be feasible, however.
Extent and fractal dimension could only identify the points where non-uniform deformation started to take place but lacked the ability to describe such a deformation. This agrees with the observation reported by Jinorose et al. 15 .

Conclusion
A novel image analysis algorithm that can be used to identify shape and describe non-uniform deformation of an input image is proposed and was tested with both computer-generated images and cubes of agar gel, which was used as the test highly shrinkable material, undergoing drying. Identifying the shape of an input image without any alignment yielded the lowest accuracy. Among the tested alignment methods, the extrema method performed adequately when both the efficiency and required computational resource were taken into consideration. The algorithm was well capable of distinguishing the different deformation patterns of the computer-generated images and, in particular, agar cubes with different solids contents. In the latter case, the developed algorithm was superior to the use of simple shape parameters, shape-parameter based algorithm and even convolutional neural networks. It is important to note that the developed algorithm can identify any shape and its change without being limited to only predetermined shapes as in the case of other existing algorithms. The developed algorithm is robust and could be further developed for real-time process control of a drying process where size and shape of a material, in particular soft material, are of concern.
Algorithm development and testing Algorithm development. Computer specification. Personal computer with Intel ® Core™ i7-10657G7 at 1.30 GHz was used in all cases. Computer is installed with NVIDIA GeForce GTX 1650 with Max-Q Design as Graphics Processing Unit with 32 GB of RAM and 512 GB of SSD. The operating system is Windows 10 Home.
Image pre-processing steps. Computer-generated images were first segmented using the algorithm described by Jinorose et al. 15 . First, artefact was eliminated from each image via the use of MATLAB ® (version R2020a, MathWork Inc., MA). Each image was then cropped to 512 × 512 pixel. Image segmentation was performed by converting an RGB image into a binary image using Otsu's thresholding method 15 . Edge detection and holes filling were subsequently performed to extract an area of interest.
As in reality a material entering an image analysis system may initially be arbitrarily oriented, a material image must first be aligned into an axis prior to being analyzed. Three alignment methods were tested in this www.nature.com/scientificreports/ study viz. equivalent ellipsis, extrema and manual alignment methods (see Fig. 1). In the case of the equivalent ellipsis method, image would be rotated until α = 0°. Note that α is an angle between the major axis of an ellipse that has the same second moment of area of interest as that of an input image and the horizontal axis 31 (Fig. 1a).
In the case of the extrema method, image would be rotated until β = 0°. In this case, β is an angle between the line connecting points p 1 and p 2 , which are the two extrema points, and the horizontal axis 28 (Fig. 1b). In the case of the manual alignment method, image would be rotated from 0° to 360° at every 1° interval. The final angle of rotation (δ) was selected as the one that gave the lowest root mean square difference (RMSD) of the distance between the centroid and the edge of the image at angle θ (d image,θ ) and that of the reference standard geometric shape (d ref , θ ) at the same angle θ (Eq. 1).
The final angle of rotation was then used to rotate the image as per the following rotation matrix (Eq. 2) using the 'affine2d' and 'imwarp' functions in MATLAB R2020a 31 : Matching input image with reference geometric shapes. Five reference geometric shapes were used in this study, i.e., triangle, rectangle, pentagon, hexagon and circle (Fig. 4). Each image itself was with the dimensions of 97 × 97 pixel, while the black background where the image sat on was with the dimensions of 512 × 512 pixel. An input image was matched to one of these geometric shapes to identify the closest shape of the image.
To identify the shape of an input image, its bounding box would be expanded until the length and width of that bounding box equal to those of the box belonging to all the reference geometric shapes. RMSD was then calculated and the most appropriate shape was selected as the one giving the lowest RMSD value.
Description of non-uniform deformation of input image. Before being able to describe non-uniform deformation of an input image, reference geometric shapes were first made to deform in various ways. The input image was first compared with the one with no deformation. The deformation patterns include swelling and shrinkage, horizontal and vertical elongations/contractions as well as convexity and concavity.
Swelling and shrinkage (or negative swelling) were achieved as per the swelling matrix (Eq. 3) 14 . Elongation, both in horizontal and vertical directions, was achieved as per the elongation matrices (Eqs. 4 and 5) 14 , while convexity and concavity functions are defined as shown in Eqs. (6) and (7) 32 . Built-in 'affine2d' and 'geomet-ricTransform2d' functions were used along with the above-mentioned equations to achieve the deformation; 'imwarp' function was used to implement all the deformation operations.
(2) Rotation matrix =  www.nature.com/scientificreports/ where S is the scaling factor. All the deformation patterns are summarized in Table 3. Parameters i, j, k, l and m are used to describe the various deformation patterns. Parameter i (Eq. 8) indicates swelling and shrinkage; '+' indicates swelling, while '−' indicates shrinkage. Parameter j (Eq. 9) indicates horizontal elongation (+) and contraction (-). Parameter k (Eq. 10) indicates vertical elongation (+) and contraction (-). Parameter l (Eq. 11) indicates convexity (+) and concavity (-). Finally, parameter m is used to indicate whether the deformation is uniform or not. This latter parameter was calculated from Eq. (13) and the aforementioned parameter l; when α is less than 0.05 and l is less than 0.01, m would indicate uniform deformation. The choices of these decision criteria is as explained in the "Results and discussion" section.
In the above equations, A is the projected area of the area of interest (AOI), L is the length of the bounding box of AOI, W is the width of the bounding box of AOI. 'major' and 'minor' are the lengths (in pixel) of the major axis and minor axis, respectively, of the ellipse that has the same normalized second central moment as that of AOI.
Comparison between developed and existing algorithms. Performance of the developed algorithm was compared with that of algorithm developed by Igathinathane et al. 17 . These investigators developed the algorithm that could identify an image into either triangle, rectangle or circle. Three shape parameters, i.e., reciprocal aspect ratio (RAR), rectangularity as well as Feret major ratio (FMR), which were calculated as per Eqs. (14-16), Table 3. Deformation patterns of reference geometric shapes. www.nature.com/scientificreports/ were used for such an identification purpose. Each of these parameters was assigned a specific range over which each of the three shapes (Shape Ig ) would be identified. Note that if the shape of an input image was far from the above three standard shapes, the algorithm would fail to capture such a shape; all peculiar shapes would be identified as an 'inclined rectangle. ' where A BB is the area of a rectangular bounding box, while D f is the maximum Feret diameter. For non-uniform deformation description, Extent (Eq. 17) 15,31 and fractal dimension (Eq. 18) 33 were also used. The description results based on these parameters were compared with those given by the newly developed algorithm.
where N r is the number of boxes intercepted with the size of each iteration of r value, while r is the size of the box.
Finally, convolutional neural networks (CNNs) with two different architectures, i.e., ResNet-50 and Xception, were tested via the use of the Deep Learning toolbox in MATLAB ® (version R2020a, MathWork Inc., MA). Stochastic gradient descent with momentum (SGDM) with InitialLearnRate of 0.01 were used 34 . The validation frequency was 315, while the MaxEpochs and MinBatchSize were both set at 10. In the cases of the presently developed algorithm and that of Igathinathane et al. 17 4200 computer-generated images were created and used to test the algorithms. On the other hand, in the case of CNNs, 45,000 computergenerated images were created and used to train (31,500 images) and validate (13,500 images) the networks. The same 4200 images were also used to test the trained CNNs.
Agar drying. Agar preparation. 2% (w/w) granulated purified agar (Product no. 1016141000, Merck Millipore Corp., Darmstadt, Germany) was mixed with either 0 or 20% (w/w) sucrose solution. The mixture was stirred at room temperature (25 ± 2 °C) at 100 rpm for 1 h, then heated to 95 °C and finally stirred at 150 rpm for 10 min. The mixture was cooled at room temperature until its temperature reached 50 °C. The cooled mixture was poured into a silicone mold to form agar cubes, each with the dimensions of 1.9 × 1.9 × 1.9 cm. The sample was allowed to set at room temperature for 1 h before being kept at 4 ± 1 °C for 24 h.
Drying experiments. Drying experiments were conducted in a hot air dryer (Memmert GmbH+ Co. KG, UM500, Schwabach, Germany) at 80 °C. Three agar cubes were taken out for moisture content determination and image acquisition at every 1 h.
Moisture content determination. A sample was weighed using a digital balance with an accuracy of ± 0.0001 g (Sartorius Lab Instruments GmbH & Co. KG, BSA224S-CW, Göttingen, Germany) and then dried in a hot air oven (Memmert GmbH+ Co. KG, UM500, Schwabach, Germany) until constant mass was obtained as per AOAC method 984. 25 (2000). The moisture content (MC) of the sample was calculated using Eq. (19): where m i is the mass of agar before drying at 105 °C, while m bd is the bone-dry mass.
Image acquisition and segmentation. The utilized image acquisition system was the one described by Jinorose et al. 15 . Only top-view images were acquired and analyzed. Image segmentation was implemented as also described by Jinorose et al. 15 .

Data availability
Experimental and relevant data are available from the authors upon suitable request. www.nature.com/scientificreports/