Loh et al. BMC Bioinformatics (2017) 18:319
DOI 10.1186/s12859-017-1720-y

S O FT W A R EOpen AccessSemi-automated quantitative Drosophila
wings measurements
Sheng Yang Michael Loh1† , Yoshitaka Ogawa2† , Sara Kawana2 , Koichiro Tamura2,3 and Hwee Kuan Lee1,3*

Abstract
Background: Drosophila melanogaster is an important organism used in many fields of biological research such as
genetics and developmental biology. Drosophila wings have been widely used to study the genetics of development,
morphometrics and evolution. Therefore there is much interest in quantifying wing structures of Drosophila.
Advancement in technology has increased the ease in which images of Drosophila can be acquired. However suchstudies have been limited by the slow and tedious process of acquiring phenotypic data.
Results: We have developed a system that automatically detects and measures key points and vein segments on a
Drosophila wing. Key points are detected by performing image transformations and template matching on Drosophila
wing images while vein segments are detected using an Active Contour algorithm.
The accuracy of our key point detection was compared against key point annotations of users. We also performed key
point detection using different training data sets of Drosophila wing images. We compared our software with anexisting automated image analysis system for Drosophila wings and showed that our system performs better than thestate of the art. Vein segments were manually measured and compared against the measurements obtained from oursystem.
Conclusion: Our system was able to detect specific key points and vein segments from Drosophila wing images withhigh accuracy.
Keywords: Drosophila, Image processing, Wing morphometrics, Automated detection

Background
Recent advances in high-throughput sequencing technology have enabled us to obtain genome informationfrom any kind of organism [1]. For humans, more thana thousand of genome sequences and single nucleotide
polymorphisms among them have been characterized [2].
This further motivated genome-wide association studies(GWAS) in the literature with the hopes of discovering thevariations responsible for genetic diseases and traits in the
genome sequences. The strategy of investigating the relationship between genotype and phenotype should succeed
in discovering the genetic basis for any inherited traits inany kinds of organisms. However, in many higher animals*Correspondence: leehk@bii.a-star.edu.sg† Equal contributors
1
Imaging Informatics Division, Bioinformatics Institute, 30 Biopolis Street,#07-01, Matrix, Singapore 138671 Singapore, Singapore
3
Research Center for Genomics and Bioinformatics, Tokyo Metropolitan
University, Hachioji, Tokyo 192-0397, Japan
Full list of author information is available at the end of the article

other than human beings, the description of phenotype
needs manpower with expertise in morphology of theanimals under consideration; such a dependency of thisapproach on available manpower may be a critical bottleneck of GWAS. Hence, development of an automatedsystem for describing phenotype may be one of the bestsolutions to this problem.
Fruit flies, in particular Drosophila melanogaster is an
important organism used for biological research fields
particularly in genetics and developmental biology. Moreover, the complete genome sequences of this species
and its related species [3] opened up an ideal opportunity for comparative genomics analyses for the systematic
understanding of phenotype and genotype relationships.
Species belonging to the genus Drosophila are particularly useful for this purpose because there are morethan 2000 described species in the genus Drosophila [4],each of which has distinct phenotypic characters. On
the other hand, this situation makes it difficult for ordinary Drosophila researchers to distinguish morphological© The Author(s). 2017 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0
International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to theCreative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Loh et al. BMC Bioinformatics (2017) 18:319

differences among different species and only a limited
number of researchers who are experts in Drosophila taxonomy can recognize key morphological characters to
identify the species.
The use of wing morphometrics for species identification have been demonstrated on a variety of insect
species [5–9], because of the simple two dimensional
structures with clear-cut patterns of veins and their crossings, which are frequently utilized as landmarks to measure distances and angles among them. There is muchinterest in quantifying wing structures of Drosophila[10–14] and Drosophila wings have been widely used to
study the genetics of development, morphometrics andevolution [15–20]. Advances in technology have increased
the ease in which images of insect wings can be acquired.
However the extraction of phenotypic information forfurther research and analysis is a manual and time consuming process. This has motivated the creation of several software which aims to automate this task. One
common method involves using software for easy annotation of landmarks on the wing (http://life.bio.sunysb.edu/
ee/rohlf/software.html, http://www.hockerley.plus.com/).
However annotating landmarks is time consuming andprone to errors [21]. More advanced software attempt
to find the location of landmarks and vein lengths using
a variety of image processing and optimization methods. These methods tend to involve fitting a set of bezier
splines onto the wing veins. WINGMACHINE [22] isa program that automatically detects and measures the
positions of veins and edges of wing blades from live flies.
MorphoJ [23] performs geometric morphometrics analysis on morphological landmark inputs. Although MorphoJ
does not detect and extract landmarks, it is a useful
tool for studying data on shape combined with molecular
genetics or ecological information. Automated systems todetect and extract information from wings from a variety
of other insect species exist [24–28]. However only a few
of these programs are able to extract specific unique key
points that are used for the prediction of the fly species.
Crnojevic et al. [28] trained a svm classifier based on Histogram of Oriented Gradient (HOG) and Complete Local
Binary Pattern (CLBP) features for vein junction detection of hoverflies. In addition the junctions were used
to construct convex hulls which were used to discriminate 4 different hoverfly species. Thus a more specific andcustomized software is needed to automate the detection
and prediction process of Drosophila wings. As mentioned previously, the wing morphology of Drosophila
species is highly favorable for automated image analysis[22]. Therefore, in this paper we introduce an automatedsystem that locates specific keypoints on a Drosophila
wing image and performs morphometrics measurementof several important vein segments. The extractionof these phenotypic data serves as a base for future

Page 2 of 14

studies and research of Drosophila wings such as specie
prediction.
The software requires the user to annotate three specific
key points on the fly wing before key points and vein segments are calculated automatically. The program provides
an intuitive user interface and users can manually edit key
points and vein segment lengths if the prediction is inaccurate. In this contribution, we report the performance ofour program and perform a comparison against WINGMACHINE [22] an existing program that automates themeasurement of Drosophila wings.

Implementation
Data set

A total of 959 flies from 16 Drosophila species wascollected in a wood on the Minami-osawa campusof Tokyo Metropolitan University in Tokyo, Japan.
The species identified according to their morphology
were Dichaetophora acutissima, Drosophila annulipes, D.
bizonata, D. busckii, D. curviceps, D. hydei, D. immigrans, D. lutescens, D. rufa, D. sternopleuralis, D. suzukii,
Hirtodrosophila sexvittata, Liodrosophila aerea, Microdrosophila sp., Scaptodrosophila coracina and Scaptomyza
graminum (Table 1). A one-sided wing was taken fromeach individual and placed in 8 μl mounting solution(ethanol : glycerol = 2 : 3) dropped on a glass slide. Aftera coverslip was placed, the slide was placed under a stereomicroscope (Nikon SMZ-10A) and the wing photograph
was taken by a digital camera (Canon EOS 40D), where theoriginal image resolutions was 1936×1288, 2816×1880,
or 3888×2592 pixels.
Table 1 Species and the number of samples used in this study
Species

FemaleMale

Dichaetophora acutissima

20

25

Drosophila annulipes

31

27

Drosophila bizonata

31

31

Drosophila busckii

35

41

Drosophila curviceps

20

26

Drosophila hydei

1120

Drosophila immigrans41

41

Drosophila lutescens

14

42

Drosophila rufa

41

41

Drosophila sternopleuralis

42

42

Drosophila suzukii

22

19

Hirtodrosophila sexvittata

41

41

Liodrosophila aerea

41

41

Microdrosophila sp.

15

10

Scaptodrosophila coracina

534

Scaptomyza graminum36

32

Loh et al. BMC Bioinformatics (2017) 18:319

A total of 600 images were used for image processing.
Figure 1 shows the arc-lengths to be measured as well as
the mean and standard deviation of the pixel intensities of
all 600 images. To conserve memory usage and improve
processing speed, small image patches around the manually annotated key points were cropped and stored insmall image patch files. In practise, only image patches
near the key point locations were needed for key point
detection. All images are also converted to gray scale
images. Let annotated images be denoted by Dm , wherem = 1 · · · N and N = 600 is the total number of images.
Key point detection

To identify a species from the wing image 13 key points
are needed. The key points are labelled ‘a’ to ‘m’ and can befound in Fig. 1. Key point ‘k’ is the boundary of the costal
fringe and the others are intersections of veins. Given a
new, unannotated fly wing image I, key points detection iscarried out in two steps. The first step aligns the annotatedimages Dm and the input image I. In the next step, key
points are detected using template matching.
To align the input image I with Dm , the user is required
to click the location of the three key points ‘a’, ‘h’ and ‘k’ on
I (see Fig. 1). A screen shot of the software can be found in
Fig. 2. An affine transformation matrix Qm is then calculated by matching the key points ‘a’, ‘h’ and ‘k’ in I and Dm .
Hence, each image in the data set Dm is associated with

Page 3 of 14

one affine transformation Qm , m = 1, · · · N. This transformation matrix is then used to transform all images inthe data set from Dm to Dm . Figure 3 shows an example of transforming Dm to Dm . In our implementation,
only pixels within a small image patch near the annotated
key points of Dm are transformed. This process greatly
increase the processing speed. As an affine transformation is used to align the images, the software is capable ofhandling fly wings of any orientation.
Template matching

Image patches between I and Dm are matched and thebest match locations are predicted to be the key point
locations. Define wI (x, y, θ) to be an image patch centred
at (x, y), oriented at an angle θ extracted from the image

I. Similarly, let wDm (x, y, θ) be an image patch extracted
from the image Dm . The objective is to find the bestmatch between wI and wDm by adjusting the image patch
parameters m, x, y, θ.
We shall describe the template matching procedureusing key point j as illustration. The procedure for matching other key points is similar. The location of key pointmj in Dm , xm
j , yj is known since Dm has been manuallymannotated. The transformed location xmcalj , yj can be


m
m
D
culated using Qm . Five templates w m xj , yj , θ areconstructed with θ = −10, −5, 0, 5, 10 degrees. Fig. 4c,

Fig. 1 Figure shows the thirteen keypoints and vein segments that can be automatically measured by the software described in this paper. In
additional, the mean and standard deviation of the pixel intensities of all 600 images used are shown as well

Loh et al. BMC Bioinformatics (2017) 18:319

Page 4 of 14

Fig. 2 Screen shot of the user interface

d, b, e and f show examples of these five templates. Next,
we want to matchthese templates to a target location

m
m
I
w xj , yj , 0 in image I. One last step before the matching process is to normalize the pixel intensities within thepatches by performing a linear transformation to make therange of pixel intensities within the image patch between
0 and 255.
m , θ
Matching is done by sliding the window w̃Dm xm
,
y
j



j

m
m
m
m
I
D
m
xj , yj , θ andwith respect to w̃ xj , yj , 0 . w̃


m
m
I
w̃ xj , yj , 0 are the normalized image patches. Figure 5
illustrates the matching process. The matching score isgiven by,
  





m
I
m
m
sm,θ (cx , cy ) = w̃Dm xmj , yj , θ − w̃ xj − cx , yj − cy , 0 (1)cx , cy is the shift of the center locations between theimage patches, with −15 ≤ cx , cy ≤ 15.  ·  is the
Euclidean norm. Finally, the predicted key point locationab

corresponds to the best match among all shifts, orientations and template images Dm .
c∗x , c∗y , θ ∗ , m∗ = arg min sm,θ (cx , cy )(2)and the key point j’s coordinate is,

  ∗

∗
xIj , yIj = xm− c∗x , ym
− c∗y
j
j

(3)

The above process is repeated for the remaining key
points on I. The details of the above algorithm canbe found in Algorithm 1 KEYPOINT-DETECTION in
Additional file 1.
Arc length calculations

Tracing the arcs as depicted in Fig. 1 is done using
active contours while the Euclidean distances between
key points are calculated trivially. In the active contour
method, an estimate of the desired curve is first generated using a template. The curve is then evolved using
a variational principle via a gradient descend method.
Two essential ingredients for applying the active contour
method are estimation of the initial curve and the existentcFig. 3 a shows an example input image I while b shows an unmapped image Dm . After applying a transformation matrix to Dm , the mapped imageDm is obtained (c)Loh et al. BMC Bioinformatics (2017) 18:319

Page 5 of 14

b

a

c

d

e

f


Fig. 4 Figure showing the target location and templates created using Dm for key point j. The target location wDm
xjm , yjm , θ = 0 (a) is created at








the location xjm , yjm on image I. c, d, b, e and f shows the 5 templates wDm xjm , yjm , θ created at xjm , yjm on image Dm with θ = −10, −5, 0,
5, 10 degrees respectivelyof a good gradient to guide the curve towards the desired
position.
Initialization of active contours

The initialization of active contour curves is critical forfinding the fly wing vein accurately. If the initial curve lies
far away from the vein, active contour might evolve thisarc to another vein.
For example, setting the initial curves to straight lines
works well for arcs fg, gh, jl and kl as these arcs are

relatively straight. This is not the case for arc ml as there
is another vein just below it. Active contour will tend to
evolve this curve towards the vein below ml.
A good initialization of curves can be obtained by mapping veins of a template fly wing image to the input image
I. This template wing should provide a good representation of wings for all different fly wing species. This
mapping is done by aligning the template fly wing image
with the input image. The alignment process has been
described in the key point detection section.

Fig. 5 The process of template matching is shown here. Each patch is moved across the search area. The Normalized Squared Summed of Pixel
Intensity Differences is calculated at each position. The position with the least pixel-wise difference corresponds to the key point location

Loh et al. BMC Bioinformatics (2017) 18:319

Figure 6 illustrates this concept. First a fly wing image
is chosen to be the template wing image (Fig. 6a). We
have selected our template wing such that key point k isequidistant to both key point j and l. Mapped arcs will
be warped badly if key point k is not equidistant to both
key points j and l. Different arcs are drawn using different
colours in Fig. 6 for clarity. Arc ml is drawn outside of thetemplate wing due to two reasons. The first is that the arc
ml of the template wing that we have chosen is quite flatand does not provide a suitable representation of other fly
wings’ ml arcs, which have more curvature. Secondly, it
can be noticed that there is another vein just below arc mlthat lies close to arc ml. Active Contour might evolve themapped arc to the wrong vein if the mapped arc happens
to lie closer to that vein and not arc ml. To prevent this,the arc ml is drawn outside the template wing image sothat the mapped arc also lies outside the fly wing. Some
results of this mapping are shown in Fig. 6b – d. One can
now notice how the mapped arc ml of Fig. 6c lies outside
the wing, rather than in between the two veins.
Generation of gradients for active contour optimization

Figure 7 shows how the input image is processed to produce a gradient for active contour optimization. The mainidea is to make the veins’ pixels very bright and for them
to be surrounded by pixels that become darker as they
become further away from the veins. A step-by-step work
flow of how the preprocessed image is obtained is shownin Fig. 7a. The original image (Fig. 7b) is blurred andinverted to obtain an image where the veins’ pixels arevery bright. The result is shown in Fig. 7c. A second image

Page 6 of 14

is needed to help produce a smooth downwards gradient away from the vein. Edges (Fig. 7d) are found on
the blurred image before applying a threshold to obtaina binary image (Fig. 7e). The binary image is dilated tomake the fly wing veins thicker before inverting the image(Fig. 7f). A distance transform (Fig. 7g) is applied before
inverting the image one last time (Fig. 7h). The two images(Fig. 7c and h) are then added pixel-wise to produce
the preprocessed image P (Fig. 7i). This process is integrated in the software and is done automatically before
performing Active Contour.
Active contour formalism

Given an image P, the arc length between 2 key points canbe found by finding the path that follows high intensity
pixels between the two points. Curves are also represented
by connecting short straight line segments. Given 2 key
points pa = (xa , ya ) and pb = (xb , yb ), n straight line segments are created to connect pa to pb by n−1 intermediate
points (x1 , y1 ). . .(xn−1 , yn−1 ) where (x0 , y0 ) = (xa , ya ) and(xn
, yn ) = (xb , yb ). The ith line segment thus has a lengthof (xi − xi−1 )2 + (yi − yi−1 )2 . An objective function for
the active contour can then be defined as,
L(x1 , y1 , . . .xn−1 , yn−1 ) = f (x0 , . . .yn , P)+αn(li −l0 )2 (4)i=1where l0 = pa − pb /n is a constant and α ≥ 0 is atuning parameter. l0 is the equilibrium length of the line
segments, i.e. if li = l0 then its contribution to the second term in Eq. 4 is zero. The function f is chosen such

a

b

c

d

Fig. 6 a shows a wing from a male Drosophila rufa fly which is being used as a template wing. This image is chosen as its key point k is equidistant toboth key point j and l. b - d shows the mapped arc lengths using the template wing for 3 different images. The colours red, blue, orange, green andmagenta represents the arc segments fg, gh, jl, kl, lm respectively

Loh et al. BMC Bioinformatics (2017) 18:319

Page 7 of 14

b

c

d

e

f

a

g

h

i

Fig. 7 Figure shows the preprocessing of image I (b) before the arc lengths are found. c and h shows the inverted image of I and the inverted
distance map of I respectively. Both images are added pixel-wise to get P, the processed image (i)that f is small along paths of high image intensity. To makethe active contours invariant with respect to fluctuationsof the overall image intensity, we make f to be invariant to image intensity multiplication, i.e. f (x0 , . . .yn , P) =f (x0 , . . .yn , aP). Hence f can be defined as:


n
P 
1
1
f (x0 , . . .yn , P) =li+
2
P(xi−1 , yi−1 ) + P(xi , yi ) + 
i=1

(5)
1
P(x0 , y0 ) + P(xn , yn )(6)2
P is the normalization on the image intensity such thatf becomes invariant under image intensity multiplication.
0 < 
1 is a small regularizer to prevent numerical
overflow. The desired curve can be found by minimizing
the objective function:using gradient descend method,

 



(t+1)
(t) (t)(t)(t)
x(t+1)
=
x
−
η∂
,
y
,
y
L
x
,
.
.
.y
k
1
n−1 (8)kkkkFinally, the arc length can then be calculated by:la,b =nx∗i − x∗i−12+ y∗i − y∗i−1

2

(9)

i=1

The details of the algorithm can be found in Algorithm
2 ARC-LENGTH in Additional file 1.

P =Species identification from fly wing images

L x∗1 , y∗1 , . . .x∗n−1 , y∗n−1 = min L(x1 , . . .yn−1 )The Drosophila species can be identified from the wing
image using 13 key points (Fig. 1). The centroid of wing
was estimated as the centroid of the octagon whose vertices were the key points ‘a’, ‘h’, ‘i’, ‘j’, ‘k’, ‘l’, ‘m’ and ‘b’. The
centroid of the octagon was computed from the areaweighted average of the centroids of the five triangles split
from the octagon. Then, Euclidean distances from thewing centroid to the 13 key points were measured and

x1 ,...yn−1(7)Loh et al. BMC Bioinformatics (2017) 18:319

Page 8 of 14

normalized by dividing by the distance between ‘a’ and ‘j’(wing length).
In addition to the normalized distances from 13 key
points to the centroid, we measured three traditional wing
indices, costal-index, c3 fringe and 5x-index [29–31] to
identify the species from a wing image. The costal-index isdefined as the ratio of the length of the vein l-m to that ofthe vein j-l. The C3 fringe is the ratio of the length of thevein k-l to that of the vein j-l. The 5X-index is the ratio of
the length of the vein g-h to that of the vein f-g. Because
arc lengths of wing veins were required to compute these
indices, we used a fitting of veins to Bezier paths. As aresult, for each wing, we obtained 16 variables to construct the training dataset and measured their values for
589 flies from 15 species excluding Microdrosophila sp.,which was used as a sample for a negative control as will
be explained later.
With this size of the training dataset, however, utilization of all of the 16 variables did not show the best
performance in the species identification in our preliminary test, where the data from 588 files were used toconstruct the training dataset and the remaining datafrom one fly was used as the test data. Therefore, we
examined the conditions that may give the best result
among all the possible numbers (2-16) and possible combinations of these numbers of variables (65519 patterns in
total) and found that the best result was obtained whensix variables, the distances from the key points e, f, g
and m to the centroid, costal-index and c3-fringe, were
employed. Accordingly, we used only these six variables
for the following discriminant analysis.
The mean vector and the variance and covariancematrix of the six variables were computed for the 15
species and for male and female separately excluding
Drosophila hydei female and Scaptodrosophila coracina
female due to a small sample size. Therefore, we obtained
the vectors for 28 groups in total. Then, the squared
Mahalanobis’ distance, D2 , [32] from each wing to eachgroup was computed. Following De Maesschalck et al.
([33]), we used D2i,j = (xj − μ̂i )T Ŝi−1(xj − μ̂i ), where xj

denotes the vector of j-th wing, μ̂i and Ŝi denote the mean
vector and the variance and covariance matrix for i-thgroup, with the aid of the R package ([34]).
For each wing image from unknown sample, theD2 values were computed with μ̂i and Ŝi for the
28 groups and the group that gave the smallest D2
value was identified to be the species of the sample.
The statistical significance of the goodness of fit to theinferred group was obtained by the chi-square test withn degrees of freedom. When the smallest D2 value waslarger than the expected value for the best-fit species
at the 0.1% level, the sample’s species was identified as‘unknown’.

To examine the efficiency of species identification by thediscriminant analysis, we used 370 additional wing images
not included in the training data from 15 species including Microdrosophila sp., which were absent in the training
data and thus expected to be identified as ‘unknown’species if the method worked correctly.
Evaluation metric of key point detection

The key points of each fly wing in the data set is foundusing the methods described above. To test the accuracyof the key point detection algorithm, 15 fly wing images,
one from each species, is manually annotated by 10 users.
This give us 10 X,Y coordinates for each key point for eachof the 15 fly wing images. The covariance matrix, which
represents the spread of the 10 annotated points for asingle key point, can be obtained from the manually annotated data. These covariances matrices are then used forbenchmarking the accuracy of the automated methods.
A total of 15 covariance matrices, 1 for each annotatedimage, are obtained for a single key point. The 15 covariance matrices are summed and averaged to obtain the
final covariance matrix for a single key point. This process
is done for all key points. Table 2 shows the covariancematrix obtained for each key point. One might noticethat the covariance matrix of key points ‘a’, ‘b’, ‘c’, ‘k’ and‘m’ contains larger values. This indicates that the manual
annotations of these key points have a large spread and
it is not easy even for humans to locate these key points
consistently.

Results
Key point detection results

Some results of Key Point Detection are shown in Fig. 8.
To determine if a predicted key point is accurate, we
check if the euclidean distance (in pixels) between the predicted key point and the ground truth is less than somethreshold. The threshold used in our experiment is 2 standard deviations away from the corresponding key point’s
covariance matrix. The results of this experiment are summarized in Table 3. The mean deviation in pixels for eachpoint is also reported. Accuracy for key points ‘a’, ‘h’ and‘k’ are not reported as these points were marked by theuser. The mean pixel deviation is also slightly higher for
key points ‘b’, ‘c’, ‘l’ and m’. These key points are locatedon slightly thicker fly wing veins as compared to other key
points. A thicker vein would increase the candidate area
for a key point and this explains the larger values in pixel
deviations.
Different sized data sets were used to predict key points.
For the data set of size 100, 100 images had key points
manually annotated and used for template matching. The
remaining 500 images were used to test the accuracy of thealgorithm. Cross validations was also done for this data
set. Similarly, for the data set of size 599, 599 images had

Loh et al. BMC Bioinformatics (2017) 18:319

Page 9 of 14

Table 2 15 images were annotated by 10 users to obtain acovariance matrix per key point per image. The averagecovariance matrix for a single key point is obtained by finding themean and standard deviation of the 15 covariance matrix
a
⎛
⎝
c
⎛
⎝
e
⎛
⎝
g
⎛
⎝
i
⎛
⎝
k
⎛
⎝
m
⎛
⎝

9.06 ± 6.37

−0.127 ± 2.64

−0.127 ± 2.64

6.20 ± 3.99

12.1 ± 14.2 −2.50 ± 3.73
−2.50 ± 3.73 2.75 ± 2.07

b
⎞ ⎛
⎠ ⎝

⎞

d
⎛

⎠

⎝

1.99 ± 0.935 0.0319 ± 0.655
0.0319 ± 0.655 2.07 ± 0.927

2.35 ± 1.48 0.263 ± 0.840
0.263 ± 0.840 2.21 ± 1.62

⎠ ⎝

⎞

h
⎛

⎠

⎝

−0.141 ± 0.710−0.141 ± 0.710

1.35 ± 0.633

6.16 ± 6.01 8.91 ± 5.29

j
⎞⎛
⎠⎝

⎞

l
⎛

⎠

⎝

13.5 ± 15.0 −3.64 ± 5.56
−3.64 ± 5.56 9.09 ± 5.76

5.22 ± 4.19 5.17 ± 2.20⎞⎠⎞

2.47 ± 1.33

−0.0481 ± 0.584−0.0481 ± 0.584

1.40 ± 0.836

1.72 ± 1.17

−0.0474 ± 0.452−0.0474 ± 0.452

1.38 ± 0.652

f⎞ ⎛2.22 ± 1.32

12.3 ± 9.36 6.16 ± 6.01

16.3 ± 13.1 5.22 ± 4.19⎠⎞

2.43 ± 1.32 0.259 ± 0.980
0.259 ± 0.980 2.06 ± 1.45

2.41 ± 0.848 0.190 ± 0.691
0.190 ± 0.691 1.81 ± 0.669

5.21 ± 2.62 −1.43 ± 1.26
−1.43 ± 1.26 1.97 ± 0.746

⎠

⎞
⎠
⎞
⎠

⎞
⎠

⎞
⎠

key points manually annotated and the remaining image
was used for testing the algorithm. We expect that thedata set of size 599 will yield better results. The accuracy results are shown in Table 3. We also downsized ourimages further and repeated the above experiments. The
results of the downsized dataset can be found in Table 4.
We also compared our method against an existing fly
wing fitting software, WINGMACHINE [22]. WINGMACHINE first detects curves in a fly wing image that fit apriori template before fitting spline curves to the image.
The location of two landmarks - the humeral break andalula notch, need to be known before curve fitting can takeplace. We ran the Wings software on ten images in ourdataset. The results are shown in Fig. 8g – i and Table 5.
WINGMACHINE tries to find the fly wing veins usingedge detection methods and results were not so ideal incases where there were background artefacts. Moreover,the edge detection algorithm sometimes fail to pick upsome veins that are too faint. This can lead to a wrong prediction on the joint locations which results in a bad fit forthe spline curves.

Arc length results

Some examples of good results of Active Contour are
shown in Fig. 8. Note that even though there are dark
patches (Fig. 8d), background artefacts (Fig. 8e) or spots
on the wing (Fig. 8f), Active Contour is still able to find thefly wing vein accurately most of the time.
We found a total of 4 bad results out of 600 images inour dataset. Users are allowed to edit bad results of arc
lengths in the fly wing software. Some examples of bad
results of Active Contour are also displayed in Fig. 8. As
can been seen in Fig. 8j and k, air bubbles are present in
the background while dust is present in Fig. 8l. It can beobserved that only the location of arc lm is incorrect andthat artefacts in these images lie above arc lm. It has been
mentioned previously that the arc lm is initialized above
the fly wing. These artefacts in the background produces
large peaks in the active contour gradient image, resulting in Active Contour getting stuck in the local maxima
of the artefact. Other arcs do not suffer from this problem
as these arcs are initialized onto the foreground which has
little artefacts. Thus only artefacts that lie just above thefly wing affects the accuracy of Active Contour.
The accuracy of our Active Contour method was also
tested. Arc lengths ground truth of fifteen different fly
wings, one for each species, are annotated manually andcompared with the arc lengths found using Active Contour. The error is calculated as the area bounded by the
annotated arc and the predicted arc. The error is thendivided by the length of the annotated arc. The results aretabulated in Table 6. This is normal as a longer arc spans alarger area of the image and this increases the amount ofnoise and background artefact that might affect the Active
Contour algorithm.
Species identification results

The efficiency of the wing image analysis for species identification is shown in Table 7. Out of 370 flies examined,
346 flies were identified as the correct species (successrate was 94%), whereas the capability of identifying correct sex was poor (success rate was 56%). This suggeststhat the wing morphology is not significantly different
between female and male, whereas the between-species
difference is significantly larger than the within-species
difference among all of the 15 species examined. In addition, it is noteworthy that the present system perfectly
recognized Microdrosophila sp., which was not includedin the training data, as ‘unknown’ species. This suggeststhat the present system has a discrimination power to give
an appropriate answer even when the species of sample is
not included in the training data. Therefore, the present
system is expected to be useful to analyze wing morphology for species identification among species belonging to
Drosophila and related genera, even after more data for
more species are added.

Loh et al. BMC Bioinformatics (2017) 18:319

Page 10 of 14

a

b

c

d

e

f

g

h

i

j

k

l

Fig. 8 a - c shows some key point prediction results of our method. d -f shows arc lengths found using our method. The colours red, blue, orange,
green and magenta represents the arc segments fg, gh, jl, kl, lm respectively. g - i shows key point and arc length predicted for the same images
using the WINGMACHINE software. Red arrows are used to indicate the errors for clarity. The errors in g and i are due to incorrect prediction of joint
locations while the error in h is due to a bad fitting results of an a priori template wing. Lastly, j - l shows three examples of bad arc length
estimation results found using Active Contour due to artefacts in the background such as air bubbles and dust particles. The red arrows show theartefacts that affect the Active contour algorithm. There are a total of four images with bad arc length estimation out of our dataset of 600 images

Discussion
Our software has managed to successfully predict landmarks on up to 15 different species of fly wing images.
Incorrect landmarks detected can be corrected using thesoftware user interface with ease. In addition, our software is able to accurately calculate and measure specificvein segments which are required for fly specie detection.
Incorrect vein segments found can also be corrected using
the user interface of the software.
In contrast to WINGMACHINE software, which makes
use of high level features to optimize fit of an a priori
model of a fly wing, our software utilizes low level features
such as edges and pixel values to find land marks and veinsegments. As mentioned by WINGMACHINE, developing detection algorithms using low levels features tend toencounter problems such as varying thickness of veins,uneven lighting and background artefacts which cause
inaccuracies during prediction. However, our software
manages to overcome the majority of these difficulties andperforms well for images in our dataset.

One advantage of using WINGMACHINE over oursoftware is that only two land marks are needed for detection instead of three. However, it should be noted that
WINGMACHINE software may require pretreatment ofimage data and this may require more manpower as compared to our software. The advantage of using MorphoJis that phenotypic data can be easily read and analysed
within the software. However, MorphoJ does not have thecapabilities to detect and extract these information fromimages automatically.
There are a few drawbacks for our algorithm. The firstis the need for a large dataset. As demonstrated in the
key point detection results section, a large dataset containing fly wing images of numerous fly species would
help increase the accuracy of our algorithm. This results
in a slow detection process as each image in the dataset
has to be loaded. To reduce the time taken for the detection process, image patches centered on the fly wing’s key
points are cropped from the original fly wing image andsaved in a separate image. This separate image is loadedLoh et al. BMC Bioinformatics (2017) 18:319

Page 11 of 14

Table 3 Results obtained using images with the original resolution
Dataset = 599
Point

Dataset = 100

1 degree rotation

5 degree rotation

1 degree rotation

5 degree rotation

Within 2
stddev (%)

Within 2
stddev (%)

Within 2
stddev (%)

Within 2
stddev (%)

Mean Pixel Dev

Mean Pixel Dev

Mean Pixel Dev

Mean Pixel Dev

b

99.3

3.13 ± 2.24

99.3

3.21 ± 2.24

98.8

3.39 ± 2.36

98.9

3.38 ± 2.39

c

86.0

4.24 ± 3.84

86.0

4.35 ± 4.32

83.5

4.57 ± 4.22

82.8

4.74 ± 4.64

d

91.5

2.59 ± 2.02

90.3

2.61 ± 2.04

87.8

2.69 ± 2.36

87.7

2.69 ± 2.41

e

92.7

2.49 ± 1.19

91.5

2.52 ± 1.23

90.1

2.70 ± 2.87

89.7

2.66 ± 2.28

f

83.7

2.56 ± 1.12

82.8

2.61 ± 2.34

78.9

2.74 ± 2.40

78.9

2.78 ± 2.50

g

95.5

2.54 ± 1.72

95.8

2.53 ± 1.67

93.0

2.65 ± 2.08

92.4

2.71 ± 2.49

i

85.0

2.72 ± 1.60

83.8

2.77 ± 1.90

79.9

2.98 ± 2.17

79.4

2.95 ± 2.19

j

90.5

2.81 ± 4.25

90.7

2.86 ± 4.24

87.7

2.97 ± 4.45

87.2

3.01 ± 4.40

l

74.0

5.45 ± 12.3

73.5

6.40 ± 15.3

65.4

11.0 ± 22.5

64.6

12.1 ± 24.6

m

94.5

5.25 ± 5.48

93.3

5.66 ± 6.62

91.0

6.37 ± 7.56

89.4

6.53 ± 7.58

The results of using different a dataset of size 599 and 100 are shown as well. As shown in the table, using a larger dataset yields a more accurate result. Also note that key
points a, h, k are excluded as these key points are annotated by the users and are thus not being predicted

instead of the whole fly wing image. As the separate image
is significantly smaller in size than the original fly wing
image, loading it is much faster and thus reduces the time
taken for key point detection. Although this method has
worked well for our dataset, it is not an efficient solution
where the size of the dataset can be very large.
The second disadvantage lies with our active contour
algorithm. The veins that are currently being measured by
active contour tend to lie on the circumference of the wingor on areas where there are no other veins lying close by.
Our active contour may face problems measuring veins
which have other separate veins in the surrounding area.
In our active contour algorithm, we first initialize the veinsegment as a straight line. In cases where there are noother veins that lie near the targeted vein, the vein segment will always evolve to the targeted vein. However, if
there exist other nearby veins that run parallel to the targeted vein, active contour may choose to evolve the veinsegment towards the other veins which are beside thetarget vein. Thus an improvement of the current active
contour algorithm is required if it is needed to measuresuch veins in the future.

Conclusion
In this paper, we have demonstrated how to locate
Drosophila wing key points and obtain arc lengths of thewing veins using image processing methods. Our results
show that our method is accurate and we compared ourTable 4 Results obtained using images 4 times smaller than the original resolution
Dataset = 599

Dataset = 100

1 degree rotation

5 degree rotation

1 degree rotation

5 degree rotation

Point

Within 2
stddev (%)

Mean Pixel Dev

Within 2
stddev (%)

Mean Pixel Dev

Within 2
stddev (%)

Mean Pixel Dev

Within 2
stddev (%)

Mean Pixel Dev

b

84.8

5.55 ± 3.53

84.2

5.65 ± 3.84

85.1

5.67 ± 3.91

85.0

5.72 ± 3.74

c

52.7

7.92 ± 5.94

53.0

7.97 ± 6.00

48.9

8.47 ± 6.33

49.0

8.45 ± 6.45

d

25.0

7.03 ± 9.16

27.2

6.99 ± 9.62

24.6

7.97 ± 10.9

26.5

7.77 ± 10.6

e

59.7

5.31 ± 6.19

62.3

5.55 ± 7.23

55.7

6.61 ± 8.76

55.4

6.91 ± 9.54

f

8.33

10.4 ± 14.4

8.17

9.84 ± 13.5

9.89

12.2 ± 16.7

10.1

12.6 ± 16.8

g

36.8

8.65 ± 9.95

37.7

9.23 ± 10.7

33.5

11.7 ± 13.5

32.8

12.0 ± 13.4

i

27.7

6.25 ± 6.73

28.7

6.24 ± 6.93

27.9

6.64 ± 7.49

27.3

6.57 ± 7.33

j

36.3

4.63 ± 5.21

34.0

4.96 ± 6.34

34.2

5.35 ± 6.11

33.1

5.76 ± 7.31

l

22.8

26.7 ± 34.2

23.0

26.3 ± 34.2

17.9

30.5 ± 35.0

17.4

30.6 ± 35.4

m

48.8

16.2 ± 13.7

47.0

16.3 ± 13.4

48.1

16.1 ± 13.8

47.6

16.1 ± 13.6

The results of using different a dataset of size 599 and 100 are shown as well. As shown in the table, using a larger dataset yields a more accurate result. Also note that key
points a, h, k are excluded as these key points are annotated by the users and are thus not being predictedLoh et al. BMC Bioinformatics (2017) 18:319

Page 12 of 14

Table 5 A total of 10 fly wings images were used for prediction
and the performance and accuracy of our method andWINGMACHINE is shown in this table. WINGMACHINE performs
poorly in some cases when the wing veins or joints are faint andare not found using edge detection methods

Table 7 Efficiency of species identification using the wing image
analysis
SpeciesNo ofsamplesNo ofsuccess
for species

(%)

No of
success
for species
and sex(%)

Dichaetophoraacutissima

3

3

100

-

-

Point Accuracy of
our method(%)

Mean pixel
deviation

Accuracy of
Mean pixel
WINGMACHINE
deviation
software method
(%)

14

100

7

50

Markedby human-

027.2 ± 10.4

Drosophila
annulipes

14

a

16

89

10

56

100.

2.59 ± 1.36

60.0

14.6 ± 16.4

Drosophila
bizonata

18

b
c

100.

2.86 ± 1.97

30.0

15.8 ± 11.2

Drosophila busckii

32

31

97

19

59

d

90.0

2.52 ± 1.08

10.0

26.1 ± 27.3

10

10

100

5

50

e

90.0

2.64 ± 1.36

40.0

19.0 ± 20.3

Drosophila
curvicepsf

90.0

3.00 ± 1.07

063.7 ± 49.0

Drosophila hydei

17

15

88

6

35

g

90.0

2.97 ± 1.32

20.0

81.8 ± 53.9

Drosophila
immigrans

38

33

87

29

76

h

Markedby human-

10.0

84.9 ± 100.

Drosophila
lutescens

20

20

100

-

-

i

60.0

2.86 ± 1.84

50.0

75.8 ± 104.

Drosophila rufa

38

35

92

20

53

j

90.0

2.43 ± 1.41

50.0

78.1 ± 108.

36

95

17

45

Markedby human

-

N.A.

-

Drosophilasternopleuralis

38

k

32

84

22

58

70.0

10.3 ± 17.6

0105. ± 106.

Hirtodrosophila
sexvittata

38

l
m

100.

3.27 ± 2.69

60.0

24.9 ± 37.8

Liodrosophila
aerea

38

38

100

25

66

Scaptodrosophila
coracina

17

16

94

11

65

Scaptomyza
graminum

24

22

92

8

33

Microdrosophila sp.

25

25

100

-

-

Total

370

346

94

179

56

Table 6 Active Contour error measurement on 15 fly wings, one
for each specie. Images are selected are random. The error is
calculated as the area bounded by the annotated arc and thepredicted arc. The error is then divided by the length of theannotated arc
Fly

ARC_F_G ARC_G_H ARC_J_L ARC_K_L ARC_L_M

D.bizonata f_1

0.495

0

0.383

1.061

3.077

D. sternopleuralis
m_2

0.831

0.07

0.838

0.624

3.491

D.annulipes_m_13

0

0.111

0.12

1.44

6.658

D.buskii_f_20

0.583

0

1.602

0.658

6.776

D.coracina_m_10

1.124

0.383

1.712

1.761

3.975

D.curviceps f_2

0.263

0.012

0.734

1.645

5.795

D.hydei_m_14

0.446

0

0.562

1.521

6.948

D.immigrans f_18

0.05

0.133

0.546

0

6.229

D.lutescens_m_10

0.539

0.086

0.527

0.328

4.607

rufa_m_20

0.153

0.021

0.985

0.712

4.8

D. suzukii_m_03

0.709

0

0.202

0.14

4.254

Di.acutissima_m_17 0.954

0

0.826

1.741

3.99

H.sexbittata_f_22

0.734

0.237

1.034

2.289

4.048

L.area_m_13

1.828

0.093

0.156

0.846

2.833

S.graminum_f_17

0.637

0.019

1.19

2.476

4.499

Hyphens indicate that the sex-separated databases are not available.
Microdrosophila sp. is not included in the training data and the success is defined as
the case of identified as ’unknown’method against WINGMACHINE, another existing fly
wing prediction software.
Listed here are some possible areas for future work. The
first is to enhance our algorithm to perform key point
predictions without the help of user defined key points.
We could also improve the arc length finding method
such that it is able to find the vein accurately even withbackground artefacts. Another possible area for future
work would be to perform further research and analysis of
the data extracted to build a fly species phenotypic tree.
Lastly, our algorithm could be extended to predict morethan just fifteen species of fly or perhaps be used to predict
the species of other insects.

Availability and requirements
Project name: DrosoWing
Project home page: http://evolgen.biol.se.tmu.ac.jp/fly/wing/Loh et al. BMC Bioinformatics (2017) 18:319

Operating system: Web browsers that support HTML5
Programming language: JavaScript, Java
Requirements: NA
License: NA

Additional file
Additional file 1: Algorithms used in the paper. This supplementary
document contains a detailed description of the algorithms used in thispaper. The Key Point Detection algorithm is shown on page 1 while the Arc
Length Detection algorithm can be found on page 2. (PDF 141 kb)
Abbreviations
CLBP: Complete local binary pattern; GWAS: Genome wide association studies;
HOG: Histogram of oriented gradientAcknowledgement
None.
Funding
Grant from Tokyo Metropolitan University. Grant from The Biomedical ResearchCouncil of A*STAR (Agency for Science, Technology and Research) Singapore.
Availability of data and materials
The web application, software and images are available for download at thissite: http://evolgen.biol.se.tmu.ac.jp/fly/wing/. A supplementary document(Algorithms.pdf) accompanies this manuscript.
Authors’ contributions
SYML and YO created the software and implemented the algorithms used in
the paper. SK and KT acquired images of Drosophila wings. KT and HKL
provided suggestions, improvements and guidance throughout the durationof this project. SYML, YO, KT and HKL contributed to the writing of themanuscript. All authors have read and approved the final version of thismanuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval
Animal experiments: The study was approved by the Research Ethics andSafety Committee of Tokyo Metropolitan University. All experiments were
performed in accordance with the Committee policies.

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.
Author details
1 Imaging Informatics Division, Bioinformatics Institute, 30 Biopolis Street,#07-01, Matrix, Singapore 138671 Singapore, Singapore. 2 Department of
Biological Sciences, Tokyo Metropolitan University, Hachioji, Tokyo 192-0397,
Japan. 3 Research Center for Genomics and Bioinformatics, Tokyo Metropolitan
University, Hachioji, Tokyo 192-0397, Japan.
Received: 8 July 2016 Accepted: 9 June 2017

References
1. Metzker ML. Sequencing technologies the next generation. Nat Rev Genet.
2010;11:31–46.
2. The 1000 G enomes Project Consortium. A map of human genome
variation from population-scale sequencing. Nature. 2010;467:1061–73.
3. Clark A, et al. Evolution of genes and genomes on the Drosophila
phylogeny. Nature. 2007;450:203–18.

Page 13 of 14

4. Markow TA, O’Grady PM. Drosophila: A Guide to Species Identification andUse. London: Academic Press; 2005.
5. Silans LMN, Passerat De, et al. Wing morphometry of Phlebotomus
perniciosus (Diptera: Psychodidae): calibration of methods with alaboratory population. Ann Trop Med Parasitol. 1996;90(5):543–50.
6. Hall MJR, MacLeod N, Wardhana AH. Use of wing morphometrics toidentify populations of the Old World screwworm fly, Chrysomya bezziana(Diptera: Calliphoridae): A preliminary study of the utility of museum
specimens. Acta Tropica. 2014;138:S49–55.
7. Francoy TM, et al. Identification of Africanized honey bees through wing
morphometrics: two fast and efficient procedures. Apidologie. 2008;39(5):488–94.
8. MORPHOMETRIC, COMPARATIVE. Preliminary study of wing morphometry
in relation to tsetse population genetics. Res Action. 2005;101:133.
9. Rohlf FJ, Archie JW. A comparison of Fourier methods for the description ofwing shape in mosquitoes (Diptera: Culicidae). Syst Biol. 1984;33(3):302–17.
10. Van Cann J, et al. Wing morphometrics as a possible tool for the diagnosisof the Ceratitis fasciventris, C. anonae, C. rosa complex (Diptera,
Tephritidae). ZooKeys. 2015;540:489.
11. Klingenberg CP, Leandro RM. Distances and directions in multidimensional
shape spaces: implications for morphometric applications. Syst Biol.
2005;54(4):678–88.
12. Dickinson MH, Hannaford S, Palka J. The evolution of insect wings andtheir sensory apparatus. Brain Behav Evol. 1997;50(1):13–24.
13. Garcia-Bellido A, De Celis JF. Developmental genetics of the venation
pattern of Drosophila. Annu Rev Genet. 1992;26(1):277–304.
14. Klingenberg CP, Zaklan SD. Morphological integration between
developmental compartments in the Drosophila wing. Evolution.
2000;54(4):1273–85.
15. Cowley DE, William RA, Rutledge JJ. Quantitative genetics of Drosophilamelanogaster. I. Sexual dimorphism in genetic parameters for wing traits.
Genetics. 1986;114(2):549–66.
16. Garcia-Bellido A. Genetic control of wing disc development in, Drosophila.
Cell patterning. Vol. 29. Amsterdam: Elsevier; 1975, pp. 161–82.
17. Diaz-Benjumea FJ, Cohen SM. Interaction between dorsal and ventral
cells in the imaginal disc directs wing development in Drosophila. Cell.
1993;75(4):741–52.
18. Garcia-Bellido A, De Celis JF. Developmental genetics of the venation
pattern of Drosophila. Annu Rev Genet. 1992;26(1):277–304.
19. Stark J, et al. The evolution and development of dipteran wing veins: asystematic approach. Annu Rev Entomology. 1999;44(1):97–129.
20. Klingenberg CP, Zaklan SD. Morphological integration between
developmental compartments in the Drosophila wing. Evolution.
2000;54(4):1273–85.
21. Dedej S, Nazzi F. Two distances of forewing venation as estimates ofwing size. J Apic Res. 1994;33(1):59–61.
22. Houle D, et al. Automated measurement of Drosophila wings. BMC Evol
Biol. 2003;3(1):1.
23. Klingenberg CP. MorphoJ: an integrated software package for geometricmorphometrics. Mol Ecol Resour. 2011;11(2):353–7.
24. Zhou Y-H, Long-Bin L, James Rohlf F. Automatic description of thevenation of mosquito wings from digitized images. Syst Biol. 1985;34(3):346–58.
25. Tofilski A. DrawWing, a program for numerical description of insect wings.
J Insect Sci. 2004;4(1):17.
26. Schroder S, et al. The new key to bees: automated identification by image
analysis of wings, Pollinating bees-the Conservation Link Between
Agriculture and Nature. Brasilia: Ministry of Environment; 2002.
27. Weeks PJD, et al. Automating insect identification: exploring thelimitations of a prototype system. J Appl Entomology. 1999;123(1):1–8.
28. Crnojevic V, et al. Image processing method for automatic discriminationof hoverfly species. Math Probl Eng. 2014;2014:Article ID 986271.
29. Gibert P, et al. Comparative analysis of morphological traits among,
Drosophila melanogaster and D. simulans: genetic variability, clines andphenotypic plasticity. Drosophila melanogaster, Drosophila simulans: So
Similar, So Different. Netherlands: Springer; 2004, pp. 165–79.
30. Sturtevant AH. The classification of the genus Drosophila, withdescriptions of nine new species, Vol. 4213. Austin: University of Texas
Publication; 1942, pp. 5–51.
31. Wheeler MR, Takada H. Diptera: Drosophilidae. Insects Micronesia.
1964;14:164–242.

Loh et al. BMC Bioinformatics (2017) 18:319

Page 14 of 14

32. Mahalanobis PC. On the generalized distance in statistics. Proc Nat Inst Sci
India. 1936;2(1):49–55.
33. De Maesschalck R, et al. The Mahalanobis distance. Chemometr Intell Lab
Syst. 2000;50:1–18.
34. R Core Team. R: A Language and Environment for Statistical Computing.
Vienna: R Foundation for Statistical Computing; 2015. https://www.Rproject.org/.

Submit your next manuscript to BioMed Central
and we will help you at every step:• We accept pre-submission inquiries• Our selector tool helps you to find the most relevant journal• We provide round the clock customer support• Convenient online submission• Thorough peer review• Inclusion in PubMed and all major indexing services• Maximum visibility for your research
Submit your manuscript at
www.biomedcentral.com/submit

BioMed Central publishes under the Creative Commons Attribution License (CCAL). Under
the CCAL, authors retain copyright to the article but users are allowed to download, reprint,
distribute and /or copy articles in BioMed Central journals, as long as the original work is
properly cited.