An image based application in Matlab for automated modelling and morphological analysis of insect wings

Despite extensive research on the biomechanics of insect wings over the past years, direct mechanical measurements on sensitive wing specimens remain very challenging. This is especially true for examining delicate museum specimens. This has made the finite element method popular in studies of wing biomechanics. Considering the complexities of insect wings, developing a wing model is usually error-prone and time-consuming. Hence, numerical studies in this area have often accompanied oversimplified models. Here we address this challenge by developing a new tool for fast, precise modelling of insect wings. This application, called WingGram, uses computer vision to detect the boundaries of wings and wing cells from a 2D image. The app can be used to develop wing models that include complex venations, corrugations and camber. WingGram can extract geometric features of the wings, including dimensions of the wing domain and subdomains and the location of vein junctions. Allowing researchers to simply model wings with a variety of forms, shapes and sizes, our application can facilitate studies of insect wing morphology and biomechanics. Being an open-access resource, WingGram has a unique application to expand how scientists, educators, and industry professionals analyse insect wings and similar shell structures in other fields, such as aerospace.

www.nature.com/scientificreports/ To the best of our knowledge, the only existing semi-automated tool for finite element modelling insect wings is our previous app, WingMesh 35 . However, WingMesh represents a significant drawback: the absence of control on the type and size of the produced mesh.
To overcome this methodological gap, here we present a powerful app, WingGram, for extracting geometric properties of insect wings and developing precise finite element wing models. We use examples of ten insect wings with different levels of geometric complexity to evaluate the applicability of WingGram (Fig. 1).

Methods and concepts
WingGram is an open-access app with a user-friendly graphical user interface (GUI) that we created in App Designer in Matlab 2020b. The fast marching algorithm 36,37 , a well-known method in computer vision, is used to extract the outer edges of an insect wing (here referred to as domain) and those of wing cells (here referred to as sub-domains) using only an input 2D image ("Extracting the boundary of the wing and wing cells"). The Ramer-Douglas-Peucker line simplification algorithm is employed to reduce the saturation of points on the detected boundaries for increasing the runtime ("Ramer-Douglas-Peucker line simplification algorithm") 38,39 . We improved our previous method, WingMesh 35,40 , and included it into WingGram to add higher flexibility in modelling 3D corrugations ("Corrugation assignment"). Python scripting in Abaqus is used to generate finite element models ("Generating models using Python scripting in Abaqus"). The fast marching algorithm ("Measurement of geometric properties") extracts the wing cell's geometric properties. The Zhang-Suen thinning algorithm 41 is used to skeletonise the lines in an input wing image and detect the location of vein junctions ("Identifying vein junctions"). The finite element simulation from Combes and Daniel 2,4,9 is used to validate generated model by WingGram ("Validation").
Extracting the boundary of the wing and wing cells. The fast marching algorithm is a recursive method in computer vision and detects an arbitrary domain 36,37 . In WingGram, we use this algorithm to detect the boundary of insect wings, their wing cells, and any discontinuity within them. Figure 2a illustrates how the fast marching algorithm works. Each square and numbers inside each figure represent a pixel and the iteration phases, respectively. Pixel 1 is an arbitrary white pixel inside the domain. When the first white pixel is selected, the code searches for white pixels around that pixel in the four cardinal directions (i.e. right, left, up and down) (Fig. 2b). When the code finds white pixels around pixel 1, it searches for the neighbour black pixels. However, this time, it checks all eight pixels around pixel 1 (i.e. Fig. 2c, pixels 2 to 9). This strategy distinguishes the neighbour cells, which might not be possible otherwise (Fig. 2d,e). The code stores the coordinates of the found white pixels and uses them as initial pixels for the next iteration. The app also keeps the coordinates of the located black pixels as the domain's boundary. In each iteration, the algorithm changes the colours of the found white and black pixels to light grey and dark grey, respectively, to avoid duplication. A similar procedure continues until there are no undetected white pixels inside the wing. Wings from ten representative insect orders were analysed and modelled using WingGram (wing drawings were taken from 16  www.nature.com/scientificreports/ Ramer-Douglas-Peucker line simplification algorithm. Using the fast marching algorithm, Wing-Gram could detect the pixels located on the boundary of a domain. However, the presence of all these pixels is not required. The higher the number of pixels at the boundary, the higher the computational complexity. Ramer-Douglas-Peucker line simplification is a recursive method that keeps only the critical points on the boundary and removes the rest 38,39 . This algorithm creates a similar curve, as the original one, with fewer numbers of pixels. In this algorithm, ε is a criterion of similarity that defines the Hausdorff distance, i.e. the maximum length, between the original curve and the simplified one 42 . A smaller ε allows more points to remain on the curve. Figure 2f illustrates how this method works on an arbitrary domain. In the first step, this algorithm marks the first and the last points of the curve (the two points overlap in a closed curve). Point P1 in the first iteration is the first found point. Then the algorithm finds the farthest point from the first kept point (Fig. 2f, iteration 2, point P2). If the distance between two found points is more than ε , the algorithm marks and keeps the new point. The algorithm then subdivides the original curve by connecting the two marked points. In each part, the algorithm finds the farthest point from the line between P1 and P2 and marks them as new fixed points if their distance from the line between P1 and P2 is more than ε (Fig. 2f, iteration 3). The algorithm recursively continues Generating models using Python scripting in Abaqus. The journal file (i.e., *.jnl file) of Abaqus is a Python script that contains all information regarding the geometry of a generated model. WingGram generates 2D-planar and 3D-shell finite element models by adapting the corresponding commands of Abaqus Python scripts. The user can define two sections for all veins and all membranes or multiple sections (""Generate Model and Figures" tab: Generating a *.jnl file"). All required commands to generate any of the two types of models are extracted and embedded in WingGram. The app automatically updates the Python scripts to develop a wing model based on the information extracted from an input image. The app uses coordinates extracted from the fast marching algorithm ("Extracting the boundary of the wing and wing cells") to define veins, membranes, and discontinuities. Python commands regarding the SHELL LOFT tool are used to assign corrugations. For this purpose, at least two loft sections and one loft path are required. The application uses extracted curved sections from the secondary image ("Corrugation assignment") as loft sections, and a linear path is defined as the loft path. The app uses the boundaries extracted from the image to determine the main domain (wing outer boundary), sub-domain (wing cells) and discontinuities of the final model ("Corrugation assignment"). After importing the model into Abaqus, the user can assign material properties to veins and membranes, set boundary conditions and loading, generate a required mesh and start simulations.

Measurement of geometric properties.
After identifying the wing's geometry, WingGram uses the obtained data to quantify wing geometric parameters, including wing cell area, length, and width. The app measures the maximum distance between two points on the boundary of a wing cell as the length. The ratio of the cell area to the cell length is the width of the cell. The cell area is determined using the method described by Bourke for irregular polygons 43 . In this method, WingGram considers a closed polygon made up of lines between N vertices x i · y i ; i is an integer between 0 and N − 1.

Identifying vein junctions.
To identify the location of the vein junction, i.e. where wing veins intersect, the Zhang-Suen thinning algorithm is used to skeletonise the veins 44 . After the skeletonisation, the width of all lines shrinks to one pixel. The app considers a pixel in the skeletonised image as a vein junction if it meets any conditions illustrated in Fig. 2l. For instance, in the first condition, if P1, P2, P4, P6, and P8 are white and P3, P5, P7, and P9 are black, P1 is a vein junction.

Validation.
To validate the methodology, we developed models of the forewing of the moth Manduca sexta and used our model to simulate the mechanical response of the wing to loading. Following the combined experimental and numerical study of Combes and Daniel 4 . we assigned homogeneous Young's moduli of 1.5 × 10 8 Nm −2 and 2.1 × 10 12 Nm −2 to membranes and veins, respectively. The wing was fixed at its base, and a point force of F = 0.003N was applied to the wing tip. We measured the maximum displacement of the wing and compared that with the experimental and numerical results of the earlier study (see the Supplementary Fig. S3).

Description of the user interface
The user interface of WingGram consists of a tab bar and a display panel. The tab-bar has four tabs, including "Home" (Fig. 3a), "Assign Corrugation" (Fig. 3b), "Morphology" (Fig. 3c), and "Generate Model and Figure"  (Fig. 3d). Below the tab bar, a display panel is embedded to show the ongoing processes (see Supplementary File S1 for the installable executive file of WingGram). The user instruction and visual description of WingGram are available in the Supplementary Video S1. "Home" tab: importing an image and detecting discontinuities and sub-domains. The "Home" www.nature.com/scientificreports/ app. Figure 2m-q shows the criteria for a suitable input image. Faded and dark images or those with salt-andpepper noises, as shown in Fig. 2m,n,q, influence the performance of the fast marching algorithm in identifying the subdomains. WingGram works particularly well if the venation patterns are clear. The application is not ideal for analysing and modelling wings with dark spots or strong pigments. Figure 2p shows thin lines in the input image, which can disturb the performance of the fast marching algorithm for differentiating neighbour subdomains. Figure 2o shows an example of a suitable image. After importing the image, the image's name and size appear in specific fields below the "Import Image" push button (Fig. 3a). By importing the image, the code automatically detects the outer boundary of the wing. In the "Discontinuity" panel, the user can mark the place of discontinuities by pushing the "Pointer" button. Activating the "Auto Detection" switch allows automatic detection of sub-domains and discontinuities. Mark a domain as the discontinuity turns its colour into dark grey. The same procedure is required to detect sub-domains using the "Sub  www.nature.com/scientificreports/ Also, there are two buttons under the list. The one with a bin icon is for deleting a domain, and the other one lets the user apply the same node density for all other domains. The same procedure is available for discontinuities. The user can switch between the list of domains and discontinuities at the top of the list. "Assign Corrugation" tab: importing the secondary image to assign corrugations. Figure 3b shows the "Assign corrugation" tab bar. The user imports the secondary image using the "Secondary Image" button ("Corrugation assignment"). A spinner is embedded to change the number of sections. We embedded two knobs in our app to change the maximum height (i.e. peak) and sharpness of corrugations. To apply changes to single sections, the user can click on a section, turn the "Link Sections" off, and then use the knobs to adjust each maximum height and sharpness parameter.
"Morphology" tab: extracting the location of vein junctions and other geometric parameters. Figure 3c shows the "Morphology" tab for extracting the geometric properties of insect wings. To use this panel, the user has to set the scale (i.e., length per pixel) and the number of the histogram bins and then push the "Start Processing" button to start extracting geometric properties from the image of the wing. As described in "Measurement of geometric properties", the result appears in several figures after the processing ends. The user can switch between the radio buttons to observe each result. Switching between the radio buttons allows visualisation of the contours/histograms of the wing cells' area, length, and width. Other radio buttons show the distributions of the vein junctions, the wing centroid's location, and the skeletonised image. On the right side of each radio button, there are two buttons. One of them saves the data of the corresponding radio button. It keeps the corresponding matrix for contours, the data of histograms, the coordinate of junctions, and the coordinate of the centroid. The second button opens the corresponding figure in a separate window. • Two-/multi-section, three-dimensional shell: These models are suitable for applying out-of-plane loadings, such as uniform pressure, impact, out-of-plane point force, and displacement 45 . • Two-/multi-section, two dimensional planar: These models are suitable for applying in-plane loadings, such as shear, tensile, or compressive loading 45 . Figure 3d shows a group of available figures. Switching between radio buttons shows figures regarding the selected option. On the "Skeletonised Image" button, the user can hide the presence of junctions. On the "Extracted Boundaries" button, the user can mark the presence of the "Main Domain," "Sub-Domains," "Discontinuities," and "Junctions" by turning off corresponding checkboxes.

WingGram in application
Ten representative wing images from different insect orders that had noticeably different shapes and venation patterns were selected to show the performance of WingGram, (Fig. 1). Wings were scaled to have the same length as 100 pixels. Figure 4 illustrates the outcome of the WingGram for each wing, including the wing's FE model, the vein junctions' location, the distributions and histograms of the cell area, cell length, and cell width (see Supplementary Model S1-S10 for FE models).
We also anticipate WingGram to apply to studies of damaged wings and fossils. Figure 5a,b show a damaged dragonfly wing and the fossil of Auroradraco eos 46 . Here we used WingGram to extract the location of the vein junctions and the area, length and width of the wing cells for both wings. Also, Using WingGram, we developed a geometric multi-section model of the damaged wing and a two-section model of the damaged wing and the fossil (find FE models in the Supplementary Models S11-S13). The models are imported into the Abaqus CAE and are meshed by triangular shell elements. In the multi-section model, we assigned different material properties to the cells. The 3D model of the damaged wing using the secondary image is generated (find the FE model in the Supplementary Model S14).
Assigning corrugations to a wing model in WingGram is convenient but might not be very accurate as the user sets them. In Fig. 5c, we used a scanning 3D measurement microscope Keyence VR 3100, to show the 3D pattern on a Sympetrum dragonfly hindwing. This image shows the main and secondary image of the wing and the 3D model generated by WingGram (find the FE model in the Supplementary Model S15). Also, two other secondary images are used for this model to show how the user can add curvature to the whole model, which is available in Supplementary Fig. S2.
WingGram is a modelling app, and virtual models generated by it could be 3D printed to construct physical models. Figure 5d shows the main and secondary image of the basal complex of the dragonfly forewing, the virtual model, and the 3D printed model. Due to the limited print area, we isolated the wing model's basal part and removed the thin membranes (See the Supplementary Model S16 for the G-Code). The isolated part of the wing was fabricated using a Prusa i3 MK3S (Prusa Research s.r.o., Prague, Czech Republic) with white coloured polylactic acid (PLA) filament (Prusa Research s.r.o., Prague, Czech Republic).
WingGram is also applicable to many 2D natural and artificial structures, such as leaves, spider nets, geographic maps and/or industrial plates and shells (see the Supplementary Fig. S1).
Validation result. Figure S3 shows the result of the FE simulations. The maximum bending displacement of our model is 5.47 mm in the z-direction, which is about 4% different from that obtained by Combes  www.nature.com/scientificreports/ www.nature.com/scientificreports/  www.nature.com/scientificreports/ The generated model in Abaqus, with its assigned boundary conditions, material properties, and loading, is available in Supplementary Model S17.

Discussion: limitations of WingGram and its advantages over other apps
WingGram equips the user with a combination of tools that can also be found in various other apps. We have also included additional unique options to WingGram, such as semi-automated finite element modelling. Here, we compare WingGram with some of the existing apps.
WingMesh. WingMesh is a Matlab-based app developed by the authors with a user-friendly interface for finite element modelling of insect wings 35  DrawWing. DrawWing extracts the location of the vein junctions in an insect wing 14 . This app needs a highresolution image (more than 2400 dpi × 2400 dpi). As mentioned by its developers, DrawWing only works with honeybee wings (Apis mellifera). In contrast, WingGram does not have any limitations regarding image resolution and is applicable for almost all types of insect wings.
NEFI & NET. These apps extract network data from an image using image processing and computer vision.
They are applicable for leaf venation, spider webs, crack paths, insect wing venations, and similar structures. They extract the location of junctions and the connection between them, which is applicable for detecting the location of vein junctions of an insect wing 18,19 . This is the only common function between the NEFI, NET, and WingGram. Our application works as accurate as NEFI and NET.

FijiWings.
FijiWings uses the advantages of the ImageJ Fiji to measure some geometric features of the Drosophila wing, like the area of the wing and trichome density. This app needs a high-resolution image to work appropriately 17 . In contrast to FijiWings, WingGram extracts the area, length, and width of subdomains, even if there are damaged input wings, and is applicable for almost all insect wings. WingGram contains some unique features applicable to field studies on insect wings regarding the wing's mechanical behaviour or morphological properties. To the best of our knowledge, WingGram is the only existing tool that can be used to develop a precise geometric model (not only mesh) using only an image of insect wings. WingGram can facilitate future studies on the insect wings and their components, such as ambient-, cross-, and longitudinal veins, corrugations, and membranes. For example, the location of junctions can be used to test fluctuating asymmetry between left and right wings, which is important for ecological studies to monitor environmental pollution 49,50 . Also, the wing shape and vein junctions' location can be used to identify insect species 14,51,52 . Extracting the area, length and width of membranes is another feature of WingGram that can enable us to study the relationship between wing form and design.
Although compared with many existing modelling methods, WingGram offers better efficiency and lesser computational time; it represents some limitations that can still be improved. WingGram models wings with homogenous thickness, while the thickness of different parts of the wing can vary from one place to another. WingGram does not model the cross-section of veins as they are in reality. The authors are currently working on developing a new version of WingGram, which can overcome these limitations and add more features to it.