Using a multiscale image processing method to characterize the periodic growth patterns on scallop shells

Abstract The fine periodic growth patterns on shell surfaces have been widely used for studies in the ecology and evolution of scallops. Modern X‐ray CT scanners and digital cameras can provide high‐resolution image data that contain abundant information such as the shell formation rate, ontogenetic age, and life span of shellfish organisms. We introduced a novel multiscale image processing method based on matched filters with Gaussian kernels and partial differential equation (PDE) multiscale hierarchical decomposition to segment the small tubular and periodic structures in scallop shell images. The periodic patterns of structures (consisting of bifurcation points, crossover points of the rings and ribs, and the connected lines) could be found by our Space‐based Depth‐First Search (SDFS) algorithm. We created a MATLAB package to implement our method of periodic pattern extraction and pattern matching on the CT and digital scallop images available in this study. The results confirmed the hypothesis that the shell cyclic structure patterns encompass genetically specific information that can be used as an effective invariable biomarker for biological individual recognition. The package is available with a quick‐start guide and includes three examples: http://mgb.ouc.edu.cn/novegene/html/code.php.

are varied, resulting in a periodic appearance of growth patterns on shells with different densities in the internal structure and chemical elements (Michio & Hiromichi, 2013). Using a combination of markrecapture methods, the features of these rings in mollusk shells can be adapted and applied to studies on the population biology of the scallop (Allison & Brand, 1995;Berkman, 1990), which are used to determine the habitat value of the ecosystems for fishery restoration and for enhancement through stocking.
Perhaps the most important tool of population biology is the ability to recognize and track individual animals over space and time.
Traditionally, this recognition has been accomplished by capturing animals and placing visible and unique marks on them (Williams, Nichols, & Conroy, 2002). Recently, molecular genetic markers such as RFLP (restriction fragment length polymorphism), RAPD (random amplified polymorphism DNA), SSR (simple sequence repeat), and SNP (single nucleotide polymorphism) have also been widely used to study the population and individual recognition (Reed, Tollit, Thompson, & Amos, 1997;Wang, 2015). However, these methods are not suitable to a larger population because of inconsistency, inconvenience and higher cost, among others. Currently, photographic mark-recapture (PMR) has gained popularity because of the advances in digital photography and image processing software. The abundance of species with variable natural marking patterns makes this an attractive method for many researchers. PMR has been employed particularly in the studies of populations of marine mammals and mammalian terrestrial predators (Fearnbach, Durban, Parsons, & Claridge, 2012;Forcada & Aguilar, 2000;Karanth & Nichols, 1998;Langtimm et al., 2004). Some image analysis algorithms have been used to extract, store, and compare pattern information from digital images (Bolger, Morrison, Vance, Lee, & Farid, 2012).
In this study, the growth patterns on the shells of Yesso scallop (Patinopecten yessoensis), Weathervane scallop (Patinopecten caurinus), and Zhikong scallop (Chlamys farreri) were analyzed using a multiscale image processing method to collect the periodic structures of growth patterns on shells for scallop individual recognition. First, we used a Gaussian-matched filter to fit the growth patterns into enhanced lines.
Next, we applied a total variation model (Rudin, Osher, & Fatemi, 1992), which was particularly effective with the enhanced image for line localization. The multiscale line maps generated by multiscale hierarchical decomposition contained a series of segmentation results at varying image resolutions of shell patterns details at different levels. This process performed an iterative segmentation at an increasing image resolution in each step, and thus detected much smaller growth rings and ribs. Finally, the periodic structures, consisting of bifurcation points, crossover points of the rings and ribs, and the connected lines, could be acquired by our Space-based Depth-First Search (SDFS) algorithm.
Meanwhile, we created a MATLAB package to implement our method for segmentation and cyclic pattern extraction of the growth rings on the scallop shell CT and digital images available in this study. The feasibility, effectiveness, and accuracy of our method were also assessed in identifying group sizes of the three species in a mixed population.

| Multiscale segmentation
The segmentation of the growth rings is the most important procedure in the identification of scallop shell rings and ribs. For the rings and ribs detection, we used matched filtering with Gaussian kernel (MFGK) ker (x, y;a, b) = −exp −(a −1 (x−b)) 2 2σ 2 , and the computed MFGK response image was as follows: where Img, (x, y), ker, a, and b denoted a shell image, a two-dimensional pixel position, two-dimensional Gaussian functions (Chaudhuri, Chatterjee, Katz, Nelson, & Goldbaum, 1989), the dilation parameter (also known as scaling parameter), and the translation parameter, respectively. r θ rotated the kernel function with an angle θ, and * represented the convolution operation in variables (x and y).
To illustrate, Figure 1 shows the MFGK response M ker (x, y; a, b) in Equation (1), which was calculated with a number of discrete values of θ from 10° up to 170° with an increment of 10° and a = 8. As shown in Figure 1, the small rings and ribs in the original image were enhanced by MFGK.
(1) M ker (x, y; a, b) = max θ (r θ (ker(x, y; a, b)) * Img(x, y)), The normalized response image was defined as follows: where μ and σ were the mean and standard deviation of the enhanced MFGK image M ker (x, y; a, b).
The multiscale hierarchical decomposition of an image f was defined as follows (Wang, Ji, Lin, & Trucco, 2013). Given an initial scale parameter λ 0 and the total variation (TV) function (Rudin et al., 1992) J(f,λ) = λ||v λ || 2 L 2 + ||u λ || BV , BV stood for the homogenous bounded total variation space equipped with the norm of total variation Based on the above enhancement with MFGK and multiscale hierarchical decomposition, many line maps u k were generated at varying image resolutions, representing different levels of line details to avoid the possible failure of feature extraction caused by a single-scale segmentation.
The multiscale segmentation results were further processed by the following procedures for better feature extraction: (1) removed the connected regions with few pixels for denoising and (2) filled the line holes based on erosion and dilation method (Gonzalez & Woods, 2007).

| Multiscale cyclic structures
To get a more accurate pattern matching, multiscale cyclic structures were extracted from the multiscale segmentation results to compute similarity between two pairs of images by our method.

| Cyclic structure extraction
The extraction procedure of cyclic structures could be achieved by the following four steps:

| Skeletonization
The binary multiscale line networks segmented in the previous step were skeletonized to identify the two-dimensional centerlines of individual branches and to determine the branch point locations of the cycle structures. A skeletonization algorithm based on morphological operations (Gonzalez & Woods, 2007) was customized for our application.

| Detection of feature points
After the skeletonization step, the feature points, which consisted of the bifurcation and crossover points in the binary skeleton line networks, were searched as follows. First, the pixels where the value sum of the 8-adjacent pixels was greater than 2 in connected components were labeled; the pixels could not be bifurcation or crossover points if the sum of its 8-adjacent pixels was less than or equal to 2. Next, the center points of all the connected components were detected as candidate feature points, and they were considered vertices (nodes), while their connected lines were considered edges (links) in graph theory. Last, the feature points were determined by removing the candidate vertices whose degree was equal to 1 iteratively, as a vertex belongs to a cycle only if it had at least two connected vertices.

| Acquisition of cyclic structures
Finding the cyclic structures from the segmented network was performed by computing the minimum cycle basis from the undirected and unweighted graph (Mehlhorn & Michail, 2009). However, because the feature points have different fixed spatial positions in the image, the spatial F I G U R E 2 Segmentation results produced by multiscale hierarchical decomposition with λ 0 = 0.01 and λ i = λ 0 2 i : (a) original CT image; (b-f) segmented image at scaling parameters λ 1 , λ 2 , λ 3 , λ 4 , and λ 5 , respectively information of the feature points and the cross-product of the vectors that form these points were used to develop the Space-based Depth-First Search (SDFS) algorithm (Algorithm 1) to find the cyclic structures.

| Description and matching of cyclic structures
To make the structure more unique, we described a cyclic structure using branch lengths, branch angles, and angles between adjacent edges, and Figure 3 shows an example description of a four-point cyclic structure. The lengths were calculated based on the pixel distance, and the angles were calculated relying on the adjacent points in the radial lines or edges. Then, we normalized all the lengths and angles using Equations (4) and (5) to preserve the scaling invariant and the guarantee rotation invariant, respectively, yielding the final feature vector given in Equation (6) as an example.
The number of feature vector elements varies with the number of cyclic feature points, and the number of branch angles varies with the style of feature points; for example, a bifurcation point usually had three branch angles, while a crossover point had four branch angles.
To obtain consistent vector length, we considered the feature points of the cycle structure as all bifurcation or crossover points, making a maximum four branch angles for each point, and the feature vector of three-, four-, five-, and six-point cyclic structures became 18, 24, 30, and 36 dimensions (the missing elements were set to 0 as demonstrated in Equation (6)), respectively.
The feature matching process can search for good similarity among all structure pairs. Matching two cyclic structures with the same number of feature points could be implemented by minimizing the similarity measure: where ṽ i and ṽ j denoted the characteristic vectors of ith and jth cycle structure in two images. The term "d(.)" denoted the measured distance between the characteristic vectors.
The essential architecture of our method for fully automated segmentation and identification of shell growth rings is shown in

| Segmentation
The results of illustrative segmentation using our multiscale PDEbased method with different scale parameters are shown in Figure 5, including CT images (top and middle) and a camera image (bottom).

| Comparison with other segmentation methods
To verify the segmentation performance of our method on the network of growth rings, we compared the present experiment results with the traditional method of image segmentation, including the Canny algorithm, morphology method, and histogram thresholding method (Gonzalez & Woods, 2007). As shown in Figure 6, our method produced the best segmentation on shell individual rings and ribs and key growing points.

| Identification of radial ribs and growth rings
The processes to find radial ribs and growth rings are described as follows. First, the segmented image with a scale was divided into six sections every 30°. In each part, a few radial ribs on the segmented image were found by the Hough transform (Gonzalez & Woods, 2007) for line detection with an angle from 0° to 30°. The radial ribs on every part were extended and crossed at the original growth point as shown in Figure 7a.
Second, we use the neighborhood-suppression method (Kovács & Julesz, 1993) for peak finding to ensure that we found spatially separated circles for the detection of growth rings as shown in Figure 7b.
Lastly, the distance of radial ribs between two consecutive growth rings, as shown in Figure 7c, was used to calculate the growth rate.

| Identification of the growth patterns
Using our SDFS algorithm, the cyclic structures consisting of bifurcation points, crossover points of the growth rings and ribs, and the connected lines were identified. Figure 8 presents multiscale segmentation and identification results of two CT images and a camera image with different scale parameters. Obviously, the more exquisitely the growth rings and ribs were segmented, the more delicate cyclic structures were presented.   Figure 9c,d,f,g,i,h, respectively, using square boxes.
We tried linear, affine, and quadratic transformation and found that the affine model was enough and robust to describe the transformation. The mosaic images aligned with the affine models of four-point, five-point, and six-point cyclic structures are shown in Figure 9e,h,k, respectively.

| Individual recognition
For each species, 50 scallops of Yesso scallop, Weathervane scallop, and Zhikong scallop were randomly selected for capturing images.
Each individual scallop had pictures taken twice in a 30-day interval, and each time, 40 segmented images were captured, and thus, 80 segmented images were produced in total for an individual scallop. Next, we divided them into two halves, with one half receiving 40 images used for training and the remaining half used for testing.
Meanwhile, we randomly selected an image as a template image for the training images for every individual scallop.

| Feature extraction
The feature set should be selected such that the between-class discrimination is maximized, while the within-class discrimination is minimized.
Indeed, in order to avoid the curse of dimensionality, it is desirable for the feature set to be as small as possible. We defined the minimum value of the similarity measure as ms = min{s ij } and used ms as a feature for individual identity. Considering that there might exist some similar cyclic structures in different-individual segmented images, an overlap ratio after registration was selected as another feature for individual identity.

| Species identification
We further tested the effectiveness of the above individual recognition method in identifying the three species Yesso scallop, Weathervane scallop, and Zhikong scallop, which are generally indistinguishable from each other by using the naked eye. Fifty scallops of each species were analyzed. As shown in Figure 11, the shapes of Yesso and Weathervane scallops, which belong to the same genus, were similar in shell shape and difficult to differentiate. Although the shape of Zhikong scallop was a little different compared to the Yesso scallop and Weathervane scallop, it was still a challenge to identify it in a big population of the three mixed species by manual sorting. Thus, we used the landmark-based geometric morphometric method (Adams, Rohlf, & Slice, 2004;Bookstein, 1993) to quantify overall shell shape for precise sorting. An important advantage of this approach is that shape information from homologous anatomical structures (landmarks), and points along curves and points on anatomical surfaces could be included in the same analysis and were termed semilandmarks. In this paper, we first obtained high-resolution bifurcation or crossover points of the shell image of each individual using our cyclic structures identification method. We then digitized the locations of five homologous landmarks on each segmented image, which are shown as red points in Figure 12a. Next, we digitized 15 semilandmarks along the ventral edge of the shell. Finally, we quantified the general shell surface by digitizing 100 cyclic structures as semilandmarks on the surface texture of each image. To accomplish this, we selected 100 cyclic structures from relatively evenly spaced surface on

Scale parameters
Classification accuracy (%)  Once all individuals were digitized, we aligned them using a generalized Procrustes superimposition (Rohlf, 1965). From the aligned shells, a set of Procrustes shape coordinates were obtained. For the pair of species, we calculated the difference in average shell shape as the Euclidean distance between species means using Procrustes shape coordinates.
We also performed a principal components analysis (PCA) to visualize the patterns of variation within and among the species. As shown in Figure 12b, the first principal component axis (PC1) and the second principal component axis (PC2) described 67.1% and 11.4% of the total variation, respectively. The two species between the Yesso scallop and Zhikong scallop or between the Weathervane scallop and Zhikong scallop could be sorted on the basis of shell shape by PC1 alone; however, there was also significant sorting between species when viewed by both PC1 and PC2. Importantly, the same genus species Yesso scallop and Weathervane scallop overlapped considerably in the morphospace, but 96% of individuals could be correctly sorted into their corresponding species using our method.

| DISCUSSION
The approach described in this paper using image processing analytical methods, which are widely used in studies on ecology and evolution (Bolger et al., 2012) This individual recognition method can be easily applied to species identification. Species identification is necessary because Yesso scallop, Weathervane scallop, and Zhikong scallop, as major economic aquaculture species, form mixed species in many areas around the China. Yesso scallop and Weathervane scallop, as cold-water species, prefer food-rich environment and have relatively poor ability to resist pollution (Shumway & Parsons, 2011). In contrast, the Asian-Pacific species, Zhikong scallop, performs better in warm and low-food availability water (Guo & Luo, 2006). One of the most compelling patterns observed in evolutionary biology is that of morphological and behavioral convergence among species inhabiting similar environments. Using our method, the Yesso scallop and Weathervane scallop, which inhabit similar environments and were more similar in their shell shape, could be easily sorted, which will be useful for scallop machine sorting in the future. recognition (Mao et al., 2013). However, our method can achieve a high-throughput operation with aid of a CT machine, an ordinary digital camera, and even mobile phones and can reduce the workload to just less than 3 hrs. Therefore, we would propose that the periodic structures of scallop could be adapted to photographic mark-recapture (PMR) to study a large population of scallops.
Genetic breeding for a higher growth rate is one of the main focuses of the scallop farming industry (Gjedrem, 1983). There are many ways to determine the growth rates in mollusk aquaculture (Michellej & Amandae, 2008). Our method of scallop individual recognition could be used to calculate the growth rate easily by characterizing shell images during a time of scallop growth and development, thus allowing the subsequent study of these growth rates as they are influenced by environmental factors or aquaculture techniques. Indeed, the establishment of a method with a possible function in predicting the growth rate in relation to scallop developmental traits could provide useful information for targeting genetic improvement of this species.

ACKNOWLEDGMENTS
This study was funded by the National Natural Science Foundation 2015ASTP-ES02).

AUTHOR CONTRIBUTIONS
Z.B., Y.W., and Z.C. conceived and designed the study. Q.X and T.W were involved in preparation of shells for image processing and conducted the major part of MATLAB package for data analysis. Z.B., Y.W., Z.C., Y.L., S.W., and L.Z. drafted the manuscript. All authors read and approved the final manuscript.

ETHICAL APPROVAL
All applicable international, national, and/or institutional guidelines for the care and use of animals were followed.
if RightHandPoints is not empty then

18.
for each point p 2m in RightHandPoints do
add p 2m determined by min(θ p 0 p 1 p 2 m ) to CyclePath

22.
elseif LeftHandPoints is not empty then

23.
for each point p 2n in LeftHandPoints do