A Study on High-Speed-Vision-BasedVibration Spectrum Imaging（高速ビジョンに基づく振動スペクトルイメージングに関する研究）

byKohei Shimasaki
島崎航平

Graduate School of Engineering
Hiroshima University
September, 2020

Contents
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.1

Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Purpose of the research . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

Outline of thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2. Related works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.1

Vibration measurement . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2

High-speed vision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92.3

Vibration source localization . . . . . . . . . . . . . . . . . . . . . . . .

10

3. Real-time high-speed vision-based vibration spectrum imaging . . . . . . . . .

13

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

3.2

Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.3

Proposed algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.4

Implementation for GPU-based high-speed vision . . . . . . . . . . . . .

17

3.5

Real-time vibration spectrum experiments . . . . . . . . . . . . . . . . .

19

3.5.1

Fixed fan with variable rotation speed . . . . . . . . . . . . . . .

19

3.5.2

Operation estimation for ﬁxed multi-copter based on the rotational speed of each propeller . . . . . . . . . . . . . . . . . . .

19

Guitar string vibration . . . . . . . . . . . . . . . . . . . . . . .

21

3.5.3
3.6

Concluding remarks. . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

4. A fast vibration source tracking algorithm using pixel-level digital ﬁlters for a
ﬂying multi-copter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

4.2

Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

4.3

Proposed algorithm of vibration source tracking algorithm using frameinterpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

i

ii

CONTENTS

4.4

The number of calculations and computational complexities . . . . . . .

33

4.5

Filter performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4.6

Execution time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

4.7

Experiments for a multi-copter . . . . . . . . . . . . . . . . . . . . . . .

46

4.7.1

Video shooting a ﬁxed multi-copter with various exposure time .

46

4.7.2

Video shooting a ﬂying multi-copter . . . . . . . . . . . . . . . .

49

4.7.3

Real-time multi-copter detection . . . . . . . . . . . . . . . . . .

51

Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

5. HFR-video-based honeybee activity sensing . . . . . . . . . . . . . . . . . . .

61

4.8

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

5.2

Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

5.3

Proposed algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

5.3.1

Pixel-level vibration source localization . . . . . . . . . . . . . .

66

5.3.2

Individual honeybee labeling . . . . . . . . . . . . . . . . . . . .

68

5.3.3

Frequency response computation for tracked honeybees . . . . .

69

5.4

Sensitivity evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

5.5

Honeybee tracking experiment . . . . . . . . . . . . . . . . . . . . . . .

75

5.6

Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

6. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
Acknowledgment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

List of Figures

1.1

Concept overview of this study. . . . . . . . . . . . . . . . . . . . . . . . .

3

3.1

Concept of pixel-level vibration spectrum imaging method . . . . . . . . . .

15

3.2

Vibration spectrum imaging for a ﬁxed fan . . . . . . . . . . . . . . . . . . . .

20

3.3

Vibration spectrum imaging for a ﬁxed multi-copter . . . . . . . . . . . . . . . .

22

3.4

Vibration spectrum imaging of guitar strings . . . . . . . . . . . . . . . . . . .

23

4.1

Concept of vibration source localization method with pixel-level digital ﬁlters 28

4.2

Flowcharts of vibration source tracking algorithms. . . . . . . . . . . . . .

30

4.3

Relationship between computational complexities and target vibration frequency (P = 3, M = 3, L = 6) . . . . . . . . . . . . . . . . . . . . . . . . .

35

4.4

Transient band in digital ﬁlter. . . . . . . . . . . . . . . . . . . . . . . . . .

37

4.5

Gain properties of elliptic ﬁlters (P = 2, 3, 5, and 10). . . . . . . . . . . . .

38

4.6

Relationship between input frame rate f0 and transition bandwidth ∆ f

. . .

38

4.7

Relationship between ﬁlter order P and transition bandwidth ∆ f . . . . . . .

39

4.8

Relationship between target frequency fN and the relative ratio of transition
bandwidth to target frequency ∆ f / fN . . . . . . . . . . . . . . . . . . . . .

39

4.9

Breakdown percentages of execution times. . . . . . . . . . . . . . . . . . .

42

4.10

Relationship between execution times and image sizes. . . . . . . . . . . .

43

4.11

Relationship between execution times and number of threads. . . . . . . . .

44

4.12

Relationship between execution times and input frame rate. . . . . . . . . .

44

4.13

Relationship between execution times and ﬁlter orders. . . . . . . . . . . .

45

iii

iv

LIST OF FIGURES

4.14

Experimental environment when video-shooting a ﬁxed multi-copter. . . . .

46

4.15

Input images when video-shooting a ﬁxed multi-copter. . . . . . . . . . . .

47

4.16

Extracted vibration features when video-shooting a ﬁxed multi-copter. . . .

48

4.17

Experimental environment for ﬂying multi-copter. . . . . . . . . . . . . . .

48

4.18

(a) Input images with tracked positions, accumulated images of (b) inputand (c) pixel-wise ﬁltered images, vibration features extracted by (d) our
proposed algorithm and (e) previous algorithm ( f0 = 2000 fps), and (e) magniﬁed images for tracked regions when video-shooting a ﬂying drone withlarge magniﬁcation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55(a) Input images with tracked positions, accumulated images of (b) inputand (c) pixel-wise ﬁltered images, vibration features extracted by (d) our
proposed algorithm and (e) previous algorithm ( f0 = 2000 fps), and (e) magniﬁed images for tracked regions when video-shooting a ﬂying drone with awider view angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

xy trajectory of extracted vibration region when a ﬂying drone was videoshot with a wider view angle. . . . . . . . . . . . . . . . . . . . . . . . . .

57

4.21

Experimental environment for ﬂying multi-copter. . . . . . . . . . . . . . .

58

4.22

(a) Input images; accumulated images of (b) input; (c) high-pass ﬁlteredimages; vibration pixels extracted by (d) our algorithm; and (e) previous
algorithm ( f0 = 2000 fps); (f) magniﬁed images. . . . . . . . . . . . . . . .

59

The xy trajectory of extracted vibration region when the video of the ﬂying
drone was shot with a wider angle. . . . . . . . . . . . . . . . . . . . . . .

60

5.1

Concept of pixel-level ﬂight activity sensing method . . . . . . . . . . . . .

80

5.2

Problems in pixel-level digital signal processing when observing fast-ﬂying
honeybees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

Input images, candidate vibration pixels, and selected ROIs for honeybees
in sensitivity evaluation. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

5.4

Magniﬁed ROI images in sensitivity evaluation. . . . . . . . . . . . . . . .

83

5.5

Averaged TFRs in the selected ROI images in sensitivity evaluation. . . . .

84

5.6

Temporal changes in image intensities in sensitivity evaluation. . . . . . . .

85

5.7

xy-trajectories and ﬂapping frequencies of labeled honeybees in 14 seconds.

85

4.19

4.20

4.23

5.3

LIST OF FIGURESv

Tracked ROI images and ﬂying activities when “bee 2” ﬂew back to its beehive for t = 6.36–6.96 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

Tracked ROI images and ﬂying activities when “bee 1” was hovering around
its beehive for t = 7.64–8.24 s. . . . . . . . . . . . . . . . . . . . . . . . .

89

Tracked ROI images and ﬂying activities when “bee 1” ﬂew away from its
beehive for t = 8.82–9.12 s. . . . . . . . . . . . . . . . . . . . . . . . . . .

91

A.1

Our prototype mobile application and future extension . . . . . . . . . . . .

99

A.2

Future onsite implementation for honeybee monitoring . . . . . . . . . . . 101

5.8
5.9
5.10

List of Tables

3.1

Execution times on a PC. . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

3.2

Execution times on PC with GPU. . . . . . . . . . . . . . . . . . . . . . .

19

4.1

Computational complexities per pixel per unit time (R = f0 / fN ) . . . . . . .

33

4.2

Execution times and image sizes (unit [s]). . . . . . . . . . . . . . . . . . .

41

5.1

Parameters used for the experiments in Subsection ??. . . . . . . . . . . . .

73

vii

Chapter 1
Introduction
1.1

Background
Vibration is a common problem in various ﬁelds such as mechanical system androtating machinery, several studies using a vibration analysis technology that grasps thetarget dynamics characteristics from actual measurement data have been conducted [1,2].
In the case of most of the exact measurement, when almost all methods are conductedusing sensors such as an accelerometer [3] and a laser Doppler vibrometer (LDV) [4],these methods are the limited number of measurement points in a narrow area or point topoint domain because of a spot measurement. These methods can’t measure the whole
dynamics characteristics of the target object, compared with vibration simulation analysis, such as ﬁnite element methods (FEM) [5]. To grasp the whole vibration dynamics bymeasuring the displacement response, vision-based methods are developed, such as the
template macthing [6, 7], optical ﬂow [8], digital image correlation (DIC) [9] and so on.
A vision sensor has potential beneﬁts as a non-contact vibration measurement in a widearea, such as numerous points measurement simultaneously. A vision sensor has highspatial resolution and distribution. The vision-based approach of displacement measurements has a limitation in terms of the spatial resolution that the image sensor of a vision
determines pixel pitch. Sub pixel-level vibration measurement with laser interferometry
using vision-based methods has been conducted such as electronic speckle pattern inter1

2

CHAPTER 1. INTRODUCTION

ferometry (ESPI) [10], as they utilize the interferometric fringe patterns for maginiﬁcationof small displacements at the laser-wavelength level. Besides, there is a limitation of alow sampling rate, which is not suitable for high-speed invisible phenomena using conventional vision systems mainly. However, they have designed standard video images forhuman’s visible ability which can see something at 30 fps (NTSC) and 25 fps (PAL) although they have a high spatial resolution. On the whole, the high-speed vision has functions of high spatial resolution and high sampling rate which can see invisible phenomena
like the acoustic signal at hundreds or thousands of fps. Because a high-frame-rate (HFR)
image sensor is regarded as a collection of photosensors that can detect temporal periodic
changes in brightness signals. We have proposed a concept of vibration features withvibration features using HFR images for capturing the time-varying changes in an image as a vibration distribution [11–13]. Real-time high-speed vision systems to capture
high-speed phenomena that the human eye cannot see have been developed at 1000 fps
level [14–17]. With the acceleration of parallel processing on ﬁeld programmable gatearrays (FPGAs) and graphics processing units (GPUs), various real-time application examples with high frame rates have been reported, such as bridge vibration analysis [18],tracking applications systems [19–21], and robot manipulation [22–25]. High-speed vision has a signiﬁcant meaning as a dynamic sensing tool which can localize and detectvibration source accurately by utilizing high spatial directivity of the image sensor, if it
can extract brightness change of extremely higher frequency level like acoustic signal
around vibration source at pixel level than the low frequency of background and so on.

1.2 Purpose of the research
In this study, we develop a vibration spectrum imaging system that implements the
pixel-level vibration frequency analysis to capture the time-varying brightness changes inimages as a vibration distribution. It realizes dynamic sensing tools for simultaneously1.2 PURPOSE OF THE RESEARCH

3

Figure 1.1: Concept overview of this study.
analyzing fast phenomena vibrating at frequencies of dozens or hundreds of Hz, whichcannot be seen by the human eye. Our high-speed vision-based vibration analysis target is a human-invisible but audible phenomenon vibrating at an audio frequency level.
The vibration spectrum imaging system consists of a high-speed vision used as a sensor
for capturing the time-varying changes in an image as a vibration distribution with theacceleration of parallel processing on GPU. We propose novel concepts to realize that avibration spectrum imaging for the high frequency targets ﬂuctuate at dozens or hundredsof Hz and high-speed moving and vibrating objects.
In our proposed vibration spectrum imaging method, we developed short-time Fouriertransform (STFT) functions for the time-varying brightness signals of all the pixels, as we
updated a conventional method [11–13] that detects only a speciﬁc frequency band witha bandpass ﬁlter and didn’t apply an unknow frequency detection using a high-frame-rate
image. These are used parallel to produce a vibration spectrum imaging that calculates
the time-frequency response in real-time for each pixel and realizes to envision the timevarying frequency response of the targets. Figure 1.1 shows the concept of our research
to visualize the time-frequency response change in the human audible phenomena.
Our concept is a pixel-level vibration spectrum imaging method that can detectsome frequencies in the audio frequency range. This concept was described in Section 3.2

4

CHAPTER 1. INTRODUCTION

in detail, implemented in the proposed algorithm of Chapter 3. To demonstrate the effectiveness of this concept, several vibration pattern recognition results are reported in
Section 3.5, with pixel-level vibration spectrum images for all the pixels of HFR images,
in which the fan, multi-copter, and guitar string are rotating or vibrating. The brightness signals for the pixels of their blades and strings ﬂuctuated at dozens or hundreds of
Hz. After we veriﬁed the eﬀectiveness of visualizing time-varying frequency responseof the targets through acoustic signal level examples such as vibrating guitar strings, wemade this concept expand two types of actual surveillance and monitoring application
for a multi-copter and a honeybee ﬂying in a cluttered outdoor scene, as we utilize apixel-level vibration source localization method.
In the case of multi-copter, a fast vibration source localization method was described in Chapter 4. Section 4.2 proposes a fast vibration source localization concept using frame-interpolated images and describes our algorithm’s details. Section 4.3 outlines
the real-time implementation of the algorithm on a high-speed vision system of 512×512
pixels at 1000 fps and evaluates the execution times of the system’s pixel-level digitalﬁlters. To verify the eﬀectiveness of our concept, Section 4.7 shows the possibility thatthe performance of the proposed algorithm is almost the same with that of the previous
algorithm through we evaluated the experiment for HFR oﬄine video and the real-time
experimental results for a ﬂying multi-copter with propellers rotating at 95 rps.
In the case of honeybee, ﬂight activity sensing was described in Chapter 5, andveriﬁed its eﬀectiveness through some experiments of ﬂying honeybees in a clutteredoutdoor scene. In Section 5.2, we propose a concept of vibration source localization forhoneybee activity sensing that executes two-step pixel-level STFTs for an HFR video toobtain the wing-ﬂapping frequencies of fast-ﬂying honeybees as their ﬂight activities, andthe details of our algorithm are described. This chapter shows the potential to quantify thewing-ﬂapping frequencies of ﬂying honeybees individually at the audio frequency level.

1.3 OUTLINE OF THESIS

1.3

5

Outline of thesis
This thesis is organized as 6 Chapters, including this introduction.
Chapter 2 summarized related works on vibration measurement, high-speed visionand vibration source tracking.
In Chapter 3, a real-time high-speed vision-based vibration spectrum imaging forvisualization of time-varying frequency response was developed to verify the eﬀectiveness of our real-time vibration spectrum imaging method. They were demonstrated using
the vibration results for a fan, multi-copter, and guitar string, which had rotating or vibrating frequencies in the range of hundreds of Hz.
Chapter 4 explains the fast vibration source localization for a multi-copter using
pixel-level digital ﬁlters was proposed to verify the eﬀectiveness of pixel-level accuracy
in vibration source tracking through several experimental results of oﬄine high-framerate videos and real-time implementation, which displayed various ﬂight scenarios of amulti-copter with its propellers rotating at 90-100 Hz.
In Chapter 5, an HFR-video-based honeybee activity sensing was proposed to realize the activity sensing of honeybees ﬂying in a natural outdoor environment by inspecting
pixel-level temporal frequency responses in the brightness in an HFR video, computed byexecuting STFTs of the brightness signals at all pixels. The eﬀectiveness of the proposed
algorithm was veriﬁed in several experiments which are how the spatial resolution of
an HFR video aﬀects the sensitivity in our honeybee activity sensing algorithm and the
trajectories and ﬂying activities of honeybees in the outdoor scene.
In Chapter 6, it summarized the contributions of this study and discussed futurework.
The appendix included the summary and future implementation of Onsite-TeamProject in the TAOYAKA program which I belonged to became a global leader with abroad perspective and learned onsite reverse innovation and bottom-up approach in a

6

CHAPTER 1. INTRODUCTION

disadvantaged area with plenty of natural resources during master and Ph.D. courses in
Hiroshima University.

Chapter 2
Related works
2.1

Vibration measurement
Vibration measurements are eﬀective, reliable, and nonintrusive technologies formonitoring and inspecting various situations such as machine operation, milling process,and structure health monitoring [26–29]. When we need to select appropriate sensors
to acquire a precise signal from a target object, we need to care about the possibility ofinterpreting what the acquired signal indicates. During the vibration transmissions, thevibration produced by a machine component consists of speciﬁc frequencies that don’tchange, although their levels vary from one location to another. Frequency analysis of
the vibration signal is widely used to diagnose machine faults [26]. In the case of cutting forces and vibration, we widely utilized these measured parameters to improve theirconﬁdence in sensor-based technologies [27]. However, it’s diﬃcult for us to detect tool
faults because the vibration dynamics of the cutting tool forces were limited. The vibration parameter is detected with a vibration sensor that can acquire an electrical signal
converted from these vibration signals by measuring a target body vibration through its
mechanical structure. We can analyze this electrical signal for vibration characteristics
such as amplitude, frequency, displacement, velocity, acceleration, phase, and period. We
can calculate a digital signal processing in the PC environment after converting the analog
signal to a digital signal.
7

8

CHAPTER 2. RELATED WORKS

Vibration measurements are commonly used such as displacement transducers (vibrometer, proximity probes), velocity transducer (velometer, LDV), accelerometer, vision
sensor.
Three types of displacement transducers are generally used for journal bearing
equipped rotating machinery for shaft vibration measurement because they preferred tomeasure the low frequency vibration less than 10 Hz. Firstly, eddy current proximity
sensor is used as a non-contact sensor to measure displacement with the electromagnetic
induction principle that changes in the resistance of the material. The eddy current proximity sensor that can measure displacement in non-contact methods is used to measure
the condition of jet engines [30]. Second, the inductive proximity switch can detect only
metal objects such as ferrous metals. When an alternating current passes through the coil,the inductance of the coil changes in present value. It is considered an optimal measurement estimated the rotor-to-stator position [31]. Third, capacitive sensors are non-contact
type sensor works on the principle of measuring the electrical capacitance of the gap
between the plates. Capacitive sensors are suited for high precision, high resolution, dynamics, and measured vibration amplitude for turbomachinery blades [32].
Velocity transducers are beneﬁtial to use high frequency range of 10 –1000 Hz,
compared with displacement transducers. LDV can especially measure the absolute velocity of a vibrating body by utilizing the Doppler shift of scattered laser light from thetarget object. As LDV utilizes the Doppler shift between the output beam and the reﬂected beam that is a frequency-modulated laser beam from the target, it can measure thetarget’ velocity. LDV is commonly used for a repeatable and constant vibration measurement such as modal analysis that is a non-destructive process of determining the inherent
dynamics properties (natural frequencies, damping ratios, and modal shapes) of vibrating
structures [4], because of speckle drop out like laser distortion of speckle noise, although
it is of high accuracy and sensitivity.
Accelerometers are electromagnetic transducers to measure vibration acceleration

2.2 HIGH-SPEED VISION

9of a target object [3]. The principle is utilizing the piezoelectric eﬀect that the force
caused by vibration causes the mass to squeeze the piezoelectric materials which produce
an electrical charge that is proportional to the force exerted upon it. Accelerometers canmeasure the vertical vibration of a target object in a various ﬁeld such as tool condition
monitoring [33,34]. Although they can measure high frequency range with a compact andlight-weight body, they need to use ampliﬁer for magniﬁcation of a weak input signal.
Vision sensor has potential beneﬁts as a non-contact vibration measurement cancapture a high spatial resolution image at the same time, such as numerous numbers of
points measurement, as above sensors are the limited number of measurement points in anarrow area because of a spot measurement. A vision sensors can measure displacement
and velocity by utilizing digital signal processing such as the template matching [6, 7],
optical ﬂow [8], digital image correlation (DIC) [9] to grasp whole vibration dynamics.
Vision sensors have a limitation of the measurable frequency range in terms of low sampling rate because input framerate is widely used at 30 fps (NTSC) and 25 fps (PAL).

2.2

High-speed vision
Conventional vision systems mainly elements standard video images for human’svisible ability which can see something at 30fps (NTSC) and 25fps (PAL) and have been
applied to various ﬁelds, such as multimedia, industry, traﬃc system, biomedical, and soon. Many appearance-based features, such as scale-invariant feature transform [35], Haarlike features [36], histogram of oriented gradients [37], and local binary patterns [38],which can be easily recognized by the human eye, have been used for target tracking invarious applications, such as traﬃc monitoring [39], human-computer interaction [40],and so on. These systems are not suitable for high-speed invisible phenomena such ascell-inspection, factory automation, and robot manipulation and so on, because of the low

10

CHAPTER 2. RELATED WORKS

framerate. Real-time high-speed vision systems to capture high-speed phenomena that thehuman eye cannot see have been developed at 1000 fps level [14–17]. We have alreadyshown the validity for high-speed tracking using high-speed vision systems such as opticalﬂow [41], camshift tracking [42], feature point tracking [43] and face tracking utilizing
Haar-like feature [44] with the acceleration of parallel processing on ﬁeld programmablegate arrays (FPGAs) and graphics processing units (GPUs). Laser triangulation [45, 46]
and LDV [47] as vibration detection using optics sensors have been proposed. Still it is
not suitable for these to measure the spatial frequency distribution simultaneously because
these approaches are mainly designated for micro-displacement measuring of one point,
such as spot inspection. On the while, in the case of high-speed vision, it can drastically
improve the amount of processing information per time. It can be visual sensing which
can give us new quality from the point of view that it can recognize and quantify highspeed phenomena human can’t see. Mainly, it has the signiﬁcant validity of a dynamic
sensing tool that can see the repetitive operation and phenomenon vibrating at soundfrequency level at dozens or hundreds of Hz as time-space distribution. And it has shown
the eﬀectiveness of cell sensing in a microchannel [48, 49] and vision-based structuralmodal analysis [50]. High-speed vision has a great meaning as a dynamic sensing tool
which can localize and track vibration source accurately by utilizing high space directivityof the image sensor, if it can extract brightness change of extremely higher frequency
level like acoustic signal around vibration source at pixel level than the low frequency ofbackground and so on.

2.3 Vibration source localization
Numerous source localization methods have been proposed [51] to simulate the human auditory sense. The binaural approach utilizes left-right-diﬀerentiated acoustic features such as interaural level diﬀerence (ILD) and interaural time diﬀerence (ITD) [52–54]

2.3 VIBRATION SOURCE LOCALIZATION

11

corresponding the human’s left and right ear. The microphone-array approach such asstereo auditory sense utilizing interaural time-intensity diﬀerence [55], correlative peakcalculation in diﬀerential arrival time distribution based on multi-microphone array system [56–58], MUltiple-Signal-Classiﬁcation (MUSIC) methods [59, 60] utilizing quadrature component in a subspace of signal and noise, and beamforming methods [61, 62]
emphasizing the directivity for designated direction, can improve localization performance and robustness by using multiple microphones organized along with a spatial geometry. They have been used in many applications such as multi-speaker recognitionin human-computer interaction [63], mobile robot audition [64], and operational vibration surveillance in product machines [65], automobiles [66], and aircraft [67]. And theopen-access software of HARK [68] that consisting of sound source localization modules, sound source separation modules, and automatic speech recognition modules has
been used widely, as they realize their actual application in computational auditory scene
analysis. Also, the many applications related to monitoring vibration [69–71] by introducing sound source localization for the manufacturing line, automobile, and plane atwork have also been reported because we want to specify vibrated abnormal parts of
those machines. Following the recent popularization of multi-copter technologies, audiosignal-based anti-multi-copter surveillance systems [72–75] can extract the sound dronepropellers emit while rotating at dozens or hundreds of Hertz, have been developed todetect illegal or abnormal drones. However, because of the low directivity of sound propagation, the localization accuracy of these methods is limited still when microphones are
distant from source objects.
Considering that a multi-copter’s propellers rotating at dozens or hundreds of Hertz
are observable in an HFR video recording at 1000 fps or more, vision-based vibration
source localization methods [12, 13] have been proposed for multi-copter and honeybee
tracking by executing digital bandpass ﬁlters for brightness signals at all pixels; they can
localize the pixels around objects vibrating at frequencies in the audio frequency range

12

CHAPTER 2. RELATED WORKS

without considering their appearances.
These vision-based methods can realize vibration source localization by detecting
the time-space change of image intensity as vibration distribution around vibration source.
This is because the image sensor is an aggregation of photodetector and image intensityof each pixel. It is time sequential signal regarding brightness change in time space andall the pixels in the image can pass the image intensity signal in a speciﬁc frequency
band by implementing pixel-level digital ﬁlter. This strong point is that it can realize thepixel-level vibration source localization only by using digital ﬁlter without implementing
the recognition algorithm based on the spatial pattern. It is a valid method in the special case that the frequency of multi copter’s propeller is much higher than one of the
backgrounds. A ﬂying multi-copter was used to demonstrate the robustness of the methods under the appearance change of occlusion, scale, focus blur, illumination condition,complex background, occlusion, scale, position, posture, and decrease resolution’s quality using low spatial resolution’s camera. Their temporal frequency responses (TFRs) inbrightness signals have been used for pixel-level vibration pattern recognition with a speciﬁc vibration pattern [76]. Since the wing-ﬂapping frequencies of ﬂying honeybees are
similar to the rotation frequencies of propellers of ﬂying multi-copters, pixel-level digital
signal processing for HFR videos could be utilized for localizing ﬂying honeybees. Theirwing-ﬂapping frequencies could be obtained as their ﬂight activities even when their appearances were sparse in the low resolution and blurred images.

Chapter 3
Real-time high-speed vision-based vibration
spectrum imaging
3.1

Introduction
Several studies that utilized vibration information for inspections in various ﬁelds

have been conducted in recent years [77–79]. Two main methods have been used to detectvibration displacement: contact measurement and non-contact measurement. In contact
measurement, contact sensors such as accelerometers [80] are directly installed on thetarget objects’ vibrated conditions. It is necessary to attach the sensor to the reference
point on the structure when attaching a contact sensor. However, it is often diﬃcult or
impossible to access the structural points at which it is installed. In non-contact measurement, the vibration displacement of the object to be measured is directly measured withsensors such as an eddy current sensor [81], an optical ﬁber sensor [82], a laser Doppler
vibrometer [83], and an oﬄine-based high-speed vision system [84]. A beneﬁcial vibration analysis of the time-varying data measured by optical sensors in the frequency domain can then be conducted after remotely measuring the vibrated conditions of the target
objects. Especially in the case of engineering, it is necessary to measure the vibration displacement [85]. In the case of a machine tool system, low vibration and noise equipment
are required for all kinds of machinery. Besides, it is always necessary to monitor and

13

14

CHAPTER 3. REAL-TIME HIGH-SPEED VISION-BASED VIBRATION SPECTRUM IMAGING

analyze the running condition of a machine during its operation [86]. It is essential toperform vibration analysis and vibration design of mechanical structures to improve theirvibration resistance. Therefore, vibration detection with non-contact sensors is more important than contact sensors when we measure and monitor vibration information simply.
Although we can see high-speed vision system as a dynamic sensing tool in subsection 2.2, various types of dynamic sensing applications have been reported using highspeed vision systems, such as the scratching behavior analysis of laboratory mice [87],cell deformation analysis using microchannels [88, 89], and vision-based structure analysis [90, 91].
Because a high-frame-rate (HFR) image sensor is regarded as a collection of photosensors that can detect temporal periodic changes in brightness signals, we have proposeda tracking algorithm with vibration features using HFR images [11]. A pixel-level vibration source localization method that uses pixel-level digital ﬁlters [12, 13] has been
proposed. A high-speed vision system is used as a sensor for capturing the time-varying
changes in an image as a vibration distribution. This method used multi-copter tracking as an example in subsection 4.1 and utilized the fact that the frequency band of themulti-copter propeller was extremely high compared to the low-frequency band, which
included the background and cyclic time series brightness signals of the propeller region.
Moreover, the presence or absence of a signal was detected using a bandpass ﬁlter witha center frequency corresponding to the rotation speed of the propeller for all the pixels.
Also, when we design to determine the target frequency band in this method, we need toknow the target frequency beforehand. Therefore, it was diﬃcult to analyze a vibrating
object when there were large frequency ﬂuctuations or a scene with an unknown vibration
frequency such as abnormal vibration.
In this vibration spectrum imaging method, we developed short-time Fourier transform (STFT) functions for the time-varying brightness signals of all the pixels, as we
renewed a conventional method [12, 13] that detects only a speciﬁc frequency band with

15

3.2 CONCEPT

input image I(x,t)

vibrated object

=

=

=

vibration
time

intensity

intensity

high-speed vision

no vibrationtime

time frequency response

vibration spectrum imaging

3RZHU

input image sequence
I(x, )I(x, )I(x, ))UHTXHQF\
3RZHU

pix
xel
el-llevel
elevelv
pixelSTFTs

peak frequency

3RZHU)UHTXHQF\

frame)UHTXHQF\

Low

High

Figure 3.1: Concept of pixel-level vibration spectrum imaging method
a bandpass ﬁlter, based on a high-frame-rate image acquired with a high-speed visionsystem. These are used in parallel to produce a vibration spectrum imaging system that
calculates the time-frequency response in real-time for each pixel and realizes to visualize
the time-varying frequency response of the targets.

3.2

Concept
The dynamic property of a vibrating object is expressed as its frequency response,

which can be obtained by converting the time-varying signal measured by an audio or
vibration sensor into a frequency in the frequency domain using the Fourier transformin the majority of sound and vibration analyses. When a vision system operates at a
suﬃciently high frame rate, temporal periodic changes can be observed in the brightness
signals at the pixels around vibrating objects corresponding to their vibration dynamic
properties in the audio frequency range. The brightness signals’ dynamic features can bequantiﬁed in the form of a pixel-level distribution by converting the brightness signals for

16

CHAPTER 3. REAL-TIME HIGH-SPEED VISION-BASED VIBRATION SPECTRUM IMAGING

all the pixels of the HFR images into their temporal frequency responses using the Fouriertransform. On the premise of a high frame rate image, the proposed concept related to
the vibration spectrum imaging based on the STFT implementation for the time-varying
signal of all the pixels is shown in Figure 4.1. When the STFT sampling rate corresponds
to the high-speed vision system’s frame rate, the STFT can be used for vibration analysis
in a frequency range of less than half the frame rate. For example, when high-speed visionat 1000 fps is used, it is possible to conduct a vibration analysis in the frequency band of
0-500 Hz. The processing ﬂow in pixel-level vibration spectrum imaging with the ﬁnal
objective of peak frequency detection is as follows.

3.3

Proposed algorithm(a) Acquisition of input HFR images
The input HFR images are acquired at time k∆t (frame number k) as I(x, y, k∆t),wherein the frame rate is f0 = 1/∆t.
(b) Computation of STFT images
For the input HFR images of K frames, the brightness signal at pixel (x, y) at timek∆t is converted to its TFR F(x, y, t) using the STFT as follows:
F(x, y, t)=(F0 (x, y, t),· · ·, F K−1 (x, y, t)),
=STFT(I(x,y,k∆t),· · ·, I(x,y,(k+K −1)∆t).

(3.1)(c) Cut-oﬀ of low-frequency components
Low-frequency components of less than fk0 = f0 k0 /(2K) are cut-oﬀ from the STFT
image F(x, y, t) as follows:
F′ (x, y, t) = (0, · · · , 0, Fk0 (x, y, t), · · · , F K−1 (x, y, t)).

(3.2)

3.4 IMPLEMENTATION FOR GPU-BASED HIGH-SPEED VISION

17(d) Normalization of cut-oﬀ STFT image
The cut-oﬀ STFT image F′ (x, y, t) is normalized as follows:′F̃ (x, y, t) = F′ (x, y, t)/|F′ (x, y, t)|(3.3)(e) Computation of peak-frequency
We detect the frequency at which the spectral component is the maximum in thefrequency band, excluding the DC component in (c) as the peak frequency for the power
spectrum calculated for each pixel.

P(x, y, t) =

3.4

f0
′
· arg max F̃k (x, y, t)
2K

(3.4)

Implementation for GPU-based high-speed vision
It is diﬃcult to realize a real-time software implementation of the above algorithm

for high-speed vision using only a CPU because the number of STFT calculations for
all the pixels becomes increasingly signiﬁcant. Therefore, in this study, we decided to
implement the real-time vibration spectrum imaging function using a GPU-based highspeed vision system that combined general-purpose computing on a graphics processing
unit (GPGPU) board and a high-speed camera that enabled parallel implementation of
the proposed algorithm. This GPU-based high-speed vision system consisted of a highspeed vision IDP Express [15], GPGPU board (Geforce GTX 1080Ti Blower, ZOTAC),
and Hewlett Packard PC (Z440 Workstation, (3.60 GHz, 6 cores, 15 MB, 2.4 GHz), Intel
C612 chipset, 32 GB memory, Windows 7 Professional OS (64 bit)). Besides, the IDP
Express consisted of the camera head and IDP Express board. This camera allowed 8 bit
gray-scale 512×512 images to be captured and processed at a maximum of 2,000 fps, andthe processed results were transferred to the PC memory at high speed. The Geforce GTX
1080 Ti (NVIDIA) on the ZOTAC Geforce GTX 1080Ti Blower had 3584 core proces-

18

CHAPTER 3. REAL-TIME HIGH-SPEED VISION-BASED VIBRATION SPECTRUM IMAGING

Image size

Table 3.1: Execution times on a PC.
128×128
256×256512×5121024×1024

K = 6443.2 ms

373.2 ms

1854.6 ms

7434.8 ms

K = 128117.7 ms

932.7 ms

4381.5 ms

17145.6 ms

K = 256

350.1 ms

2280.3 ms

9988.3 ms

40039.8 ms

K = 5121028.9 ms

5178.1 ms

21860.8 ms

84388.7 mssors operating at a basic operating frequency of 1480 MHz, as well as 11 GB of GDDR
5X global memory with a bandwidth of 484.4 GB/ s, which provided the computational
performance needed to execute ﬂoating point arithmetic with 11.3 Tﬂops. We evaluate
the execution times of our algorithm using the personal computer (PC) and GPU to calculate the subprocesses (a)–(e) in the pixel-level vibration spectrum imaging. Table 3.2
summarizes our algorithm’s execution times for various image sizes, where K indicatesthe number of frames in the STFT computation. Our algorithm’s execution time increased
in proportion to the total number of image pixels and the number of images in the STFT
computation. In the case of the real-time software execution, the operable frame ratesof the vision system were 840, 206, 56, and 13 fps for images with sizes of 128×128,
256×256, 512×512, and 1024×1024, respectively, when K = 128. The software could
process the low-resolution images in real-time at several frames per second. In contrast,
our algorithm should be accelerated for the real-time processing of the higher resolution
images at dozens of frames per second by implementing the parallel processing logicof our algorithm on speciﬁc accelerators such as general-purpose GPUs. When we implemented a pixel-level STFT calculation for K = 128 frames using 512×512 images at2000 fps, the execution could be performed in 17.89 ms. It was possible to conduct all of
the processes from the pixel-level STFT calculation to the peak frequency image output.
We conﬁrmed that the vibration spectral image could be displayed in real-time, including
a series of operations, at 50 fps or higher.

19

3.5 REAL-TIME VIBRATION SPECTRUM EXPERIMENTS

Table 3.2: Execution times on PC with GPU.
Image size 128×128
256×256512×512 1024×1024
K = 640.66 ms

2.20 ms

7.60 ms

31.06 ms

K = 128

1.19 ms

4.85 ms

17.89 ms

74.08 ms

K = 256

2.55 ms

10.92 ms

42.06 ms173.90 ms

K = 5125.69 ms

24.50 ms

95.93 ms

391.85 ms

3.5

Real-time vibration spectrum experiments

3.5.1

Fixed fan with variable rotation speed
First, we conducted a real-time vibration imaging experiment using 512×512 im-ages at 2000 fps for a ﬁxed fan whose blade rotation speed ﬂuctuated. In the center of theimage, we ﬁxed a fan with four 5 mm long two-bladed propellers. We then changed their
speed to 275, 450, 150, and 350 rotations/second at time t = 0, 5, 10, and 15 s, respectively. We considered 128 frames (K = 128) in this experiment when computing the STFT
operations;the cut-oﬀ frequency was fk0 = 30 Hz. Figure 3.2 shows the (a)input image;(b) peak-frequency image, (c) STFT results for t = 0, 5, 10, and 15 s; and (d) vibration
frequency spectrum for a 1.5 s period at 15 s. We can see peak frequency ﬂuctuated fromvibration frequency spectrum around a fan in (c) and (d). Based on these results, evenwhen the target vibration frequency ﬂuctuated, a peak frequency corresponding to twice
the rotational speed of a fan was detected, and pixel-level vibration spectrum imaging was
realized in real-time.

3.5.2

Operation estimation for ﬁxed multi-copter based on the rotational speed of each propeller
We conducted a real-time vibration imaging experiment using 512×512 images at

500 fps for a multi-copter whose propeller rotation speed ﬂuctuated. In the center of theimage, we ﬁxed a multi-copter (RC EYE One Xtreme, CEI Conrad), 225×225×80 mm

20

CHAPTER 3. REAL-TIME HIGH-SPEED VISION-BASED VIBRATION SPECTRUM IMAGING

t=0s

t=5s

t = 10 s

t = 15 s(a) input images

t=0s

t=5s

t = 10 s

t = 15 s

0.10
0.08
0.06
0.04
0.02
0.00
0

200

400

600

800

1000

0.12

900 Hz

0.10
0.08
0.06
0.04
0.02
0.00
0

200

400

600

800

1000

0.08

Power Spectrum

550 HzPower Spectrum

Power Spectrum

Power Spectrum

0.12

300 Hz

0.07
0.06
0.05
0.04
0.03
0.02
0.01
0.00
0

200

400

600

1000

30

(b) peak-frequency images

800

0.14
0.12

700 Hz

0.10
0.08
0.06
0.04
0.02
0.00

1000

0

200

400

600

800

Frequency [Hz]Frequency [Hz]Frequency [Hz]Frequency [Hz]

t=0s

t=5s

t = 10 s

t = 15 s

1000(c) STFT results (around a fan)
[Hz]
1000

0

0.1

750
500

250
0
0

0.5

1.0

1.5 [s](d) vibration frequency spectrum for a 1.5 s period at 15 s

Figure 3.2: Vibration spectrum imaging for a ﬁxed fan

in size, with four 138 mm long two-bladed propellers. At time t = 0, 5, 10, and 15 s, thepropeller operation was changed to hovering, moving forward, turning left, and turning
right. In the experiment, the STFT computation was conducted for K = 128 frames;3.6 CONCLUDING REMARKS21

the cut-oﬀ frequency was fk0 = 30 Hz. Figure 3.3 shows the (a)input image; (b) peakfrequency image for t = 0, 5, 10, and 15 s; and (c) moving direction classiﬁcation basedon the peak frequency. These results could be used to estimate the moving direction of themulti-copter from the top view, and pixel-level vibration spectrum imaging was realized
in real-time.

3.5.3

Guitar string vibration
We conducted a real-time vibration imaging experiment with a vibrating open guitar

string using 512×512 images at 1000 fps. We considered 128 frames (K = 128) in thisexperiment when computing the STFT operations; the cut-oﬀ frequency was fk0 = 0 Hz.
Figure 3.4 shows the (a)input image and (b) peak-frequency image for t = 0, 5, 10, and
15 s. These results indicate that peak frequencies of approximately 82, 110, 147, 196, 247,
and 330 Hz, corresponding to the frequencies of the 6th to 1st strings, were sequentiallydetected from the top, and real-time imaging was realized for objects with vibrations at
the audio frequency level.

3.6

Concluding remarks
In this chapter, we developed a real-time vibration spectrum imaging system inwhich STFT functions were implemented parallel for the time-varying signals of all the
pixels. This camera allowed vibration imaging at 50 fps or higher using 512×512 images
at hundreds and thousands of frames per second to be realized in real-time. Also, we
demonstrated its eﬀectiveness through veriﬁcation experiments using several vibrating
objects. We are planning to research mechanical vibration analysis, structural inspection,and biological dynamics inspection as applied research on real-time vibration spectrum
imaging.

22

CHAPTER 3. REAL-TIME HIGH-SPEED VISION-BASED VIBRATION SPECTRUM IMAGING

t=0s

t=5s

t = 10 s

t = 15 s(a) input images

t=0s

t=5s

t = 10 s(b) peak-frequency images

riset = 15 s
30

250

descent
30(c) moving direction classification based on the peak frequency
Figure 3.3: Vibration spectrum imaging for a ﬁxed multi-copter

250

23

3.6 CONCLUDING REMARKS

t=0s

t=5s

t = 10 s

t = 15 s(a) input images

t=0s

t=5s

t = 10 s(b) peak-frequency images
Figure 3.4: Vibration spectrum imaging of guitar stringst = 15 s
0350Chapter 4
A fast vibration source tracking algorithm using
pixel-level digital ﬁlters for a ﬂying multi-copter
4.1

Introduction
With the recent decrease in the price of multi-copters, the use of drone technology

has expanded globally. Drone technology has been applied in various aerial photography
applications for the press, entertainment, infrastructure maintenance, and disaster monitoring. It is predicted that a vast worldwide market for drone technology will emerge
in the next decade, including many aerial services such as delivery drones to transport
packages and wireless communication with drone-based mobile stations as well as aerial
photography. However, with the expansion of drone technology, the number of untoward
or unfortunate incidents with multi-copters caused by out-of-control crashes, terrorism,
smuggling, and privacy invasion has rapidly increased. These illegal and immoral intrusions with multi-copters have become a severe social problem when promoting drone
technology [92].
As a countermeasure against such illegal intrusions by multi-copters, there is a fastgrowing demand for surveillance technologies that can detect, localize, and identify multiple types of multi-copters ﬂying in the air and regulations for air traﬃc control andoperations involving multi-copters. Acoustic signal-based drone monitoring systems that25

26 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

can identify the sound at dozens or hundreds of Hz generated by the rotating propellerof a ﬂying multi-copter have been developed [72–75]. There have been several systems
introduced for drone monitoring at airports and infrastructures in the European Union andthe United States of America. These acoustic-based systems can detect the presence of amulti-copter by acoustic signal processing. However, they remain inaccurate in localizing and tracking a multi-copter when they are ﬂying outdoors at a long distance from theacoustic sensors because of the very limited spatial directivity of acoustic sensors.
Optical image sensors have an advantage over acoustic sensors with their acutespatial directivity compared with acoustic sensors. Moreover, many vision-based target
tracking systems that can accurately localize a target object in images have been developed for robot vision [93]. These systems used conventional image sensors with standard video signals at dozens of frames per second and were suitably designed for the
naked human eye. Many appearance-based features, such as scale-invariant feature transform [35], Haar-like features [36], histogram of oriented gradients [37], and local binary
patterns [38], which can be easily recognized by the human eye, have been used for target
tracking in various applications, such as traﬃc monitoring [39], human-computer interaction [40], and so on. To detect ﬂying objects, Rozantsev et al. [94] proposed an integrated
method that utilized both motion estimation from image sequences and appearance-based
features. It is diﬃcult to detect time-varying brightness signals at hundreds of Hertz inimage sequences with dozens of frames per second, which correspond to the rotating
propellers of a multi-copter. The appearance-based method has limitations in tracking
accuracy and robustness when the appearance of a ﬂying multi-copter changes according
to the lighting condition, background, occlusion, scale, and so on.
As we introduced a pixel-level vibration source localization method without any
knowledge of its appearance features in subsection 3.1, Jiang et al. [12, 13] proposed
this method to extract time-varying brightness changes with digital ﬁlters at all the pixels
in HFR images. However, its computational complexity was too high for the real-time4.2 CONCEPT

27

execution of digital ﬁlters at all the pixels in the high image resolution of HFR images.
They succeeded in implementing ROI image tracking in real-time, as they utilized small
resolution in a limited area, compared with input resolution because of the limitation of
execution time. Thus, we propose a novel fast vibration source tracking algorithm thatcan reduce the computational complexity of pixel-level digital ﬁlters for HFR videos byvirtually adjusting the sampling rates of the digital ﬁlters to twice the target’s vibration
frequency with frame interpolation through sampling rate locks in the object’s frequency.
Then we also developed a real-time multi-copter tracking system on which pixel-level
digital ﬁlters can detect vibration distributions around the rotating propellers of a ﬂying
multi-copter. Finally, we realized that high-speed and high spatial resolution camera could
speed up. We need to implement the digital ﬁter regarding the higher resolution to monitor
vast areas.

4.2

Concept
Assuming that the frame rate of a vision system is a much higher than the target’s

vibration frequency, and its image sensor is a collection of photodetectors, Jiang’s method[12] brought the realization that vibration source localization by extracting the speciﬁcfrequency components of brightness signals at pixels around vibrating objects with pixellevel digital ﬁlters. Figure 4.1 illustrates this concept. The advantage of this concept isthat the simple implementation of digital ﬁlters on all the pixels enables vibration source
localization without appearance-based recognition. This method is eﬀective in localizing
a vibrating object, whose vibration frequency is much larger than that of its backgroundscene, such as a ﬂying multi-copter with rotating propellers, because of its robustness totime-varying appearances that depend on lighting conditions, lens blur, scale, and posturechanges, and occlusion.
Real-time multi-copter tracking with pixel-level digital ﬁlters has been reported for28 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

input image ( , )
= 2

= 1

= 3

vibra!ng object( , )( , 2)( , 1)

input image

vibra!on!me

pixel-level
digital ﬁltersn frameintensity

intensity

high-speed vision

no vibra!on!mevibra!on-based
image features(posi!on, size,..)vibra!on distribu!onFigure 4.1: Concept of vibration source localization method with pixel-level digitalﬁltersimages at 1000 fps [13]. In contrast, the number of pixels on which digital ﬁlters were
implemented for real-time execution has been limited due to its computational complexity. Digital ﬁlters were implemented at 128×128 pixels selected from 512×512 input
images in [13]. Thus, it is crucial to accelerate the speed of a vibration source localization algorithm that simultaneously executes digital ﬁlters at much more pixels in HFR
images. This method enables real-time tracking with much higher spatial resolution of amulti-copter ﬂying in wide-view scenarios.

4.3

Proposed algorithm of vibration source tracking algorithm using frame interpolation
In general, we can design the digital ﬁlter passing a speciﬁed vibration frequency

component if the object frequency is lower than the Nyquist frequency, half of the framerate. On the while, the computational complexity of the proposed method [11–13] depends

4.3 PROPOSED ALGORITHM OF VIBRATION SOURCE TRACKING ALGORITHM USING FRAME INTERPOLATION29

on the framerate. It calculates the bandpass ﬁlter, which refers to center frequency andobject frequency by each pixel for video image recording at ﬁxed framerate. This is why
it is not necessarily eﬃcient for the previous method to calculate extra complexity in thecase of tracking objects whose frequency is much lower than the framerate.
We propose the high-speed vibration source tracking algorithm which reduced thecomplexity per unit time by implementing the pixel-level high pass ﬁlter that can pass thefrequency component around target frequency. After we set framerate twice as many astarget frequency using downsampling by frame interpolation, assuming we have alreadyknown the target frequency, which is tracking the object’s frequency. Figure 4.2 shows
the ﬂowchart of this method comparing the previous method. We selected the Linear
frame interpolation method, which can equalize the frame period, virtually not selecting
the Nearest neighbor interpolation method, which has the possibility of a frame period
being irregular. It has the purpose of accurately detecting the frequency component byusing a digital ﬁlter, assuming the sampling period is uniform. The detailed processing
ﬂow of this proposed algorithm will show from now. The vibration feature to be evaluatedin this study is calculated as follows:(1) Image acquisition
From an HFR camera operating at a ﬁxed frame rate f0 , an input image at timet = k∆t is captured as follows:

Ik (x) = I(x, k∆t),(4.1)

where k, ∆t0 = 1/ f0 , and x = (x, y) are the frame number, frame interval of the HFRcamera, and coordinate value of pixels, respectively.
(2) Downsampling with frame interpolation
When the target’s vibration frequency is fN (t) at time t, the virtual sampling time T K

30 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

frame rate

frame rate(1) Image acquisi!on0(1) Image acquisi!on
0

( )

( )

(2) Down-sampling with
Frame interpola!oñ ( )(3) Pixel-level HPF(2) Pixel-level BPF

( )

( )

(4) Oﬀset processing with
max-min diﬀerence(3) Oﬀset processing with
max-min diﬀerence

( )

( )

( )

( )

(5) Image accumula!on(4) Image accumula!on

( )

( )

( )

(6) vibra!on region extrac!on( )(5) vibra!on region extrac!on
( )

( )

(a) our proposed method(b) previous method

Figure 4.2: Flowcharts of vibration source tracking algorithms.
˜ T K ) is updated by adding twice the cycle time ofof the downsampled video I˜K (x) = I(x,the target vibration T N (T K−1 ) = 1/ fN (T K−1 ) to the previous sampling time T K−1 as follows:
T K = T K−1 + 2T N (T K−1 ),(4.2)

˜ T K ) are generated by executing linear
Downsampled images at time T K , I˜K (x) = I(x,frame interpolation between two input images Ik′ (x), Ik′ +1 (x), which are sampled at the
two nearest times around time T K as follows:˜IK (x) = (k′ +1)− T K Ik′ (x)+ T K −k′ Ik′ +1 (x),
∆t
∆t

(4.3)

where k′ is an integer that satisﬁes k′ ∆t < T K ≤ (k′ + 1)∆t.
(3) Pixel-level high-pass ﬁlters
By applying inﬁnite impulse response (IIR) high-pass ﬁlters to time-varying brightness signals at all the pixels x in the downsampled images I˜K (x), the high-pass ﬁltered4.3 PROPOSED ALGORITHM OF VIBRATION SOURCE TRACKING ALGORITHM USING FRAME INTERPOLATION31

images gK (x) = g(x, T k ) at time T K are computed as follows:gK (x) =P−1
X

b s I˜K−s (x) −

s=0

P−1
X

a s gK−s (x),

(4.4)

s=1

where P is the ﬁlter order, and a s , b s are the tap coeﬃcients; they determine the propertyof the high-pass ﬁlter. I˜K (x) and gK (x) are initially set to zero images. The symmetry
frequency component has an aliasing eﬀect on the downsampled images, concerning the
Nyquist frequency fN , enabling a high-pass ﬁlter with a pass frequency (1 − a/2) fN toequivalently function as a bandpass ﬁlter, whose center frequency and bandwidth are fN
and a fN , respectively.
(4) Oﬀset processing
To remove the oﬀset values in the downsampled images and their high-pass ﬁlteredimages, the oﬀset-diﬀerence images IKA (x) = IKA (x, T K ) and gAK (x) = gAK (x, T K ) at time T K
are computed by diﬀerentiating the maximum and minimum values of I˜K (x) and gK (x)for time t = T K − T M to T K as follows:
IKA (x) = IKM (x) − IKm (x),(4.5)

gAK (x) = gKM (x) − gmK (x).

(4.6)

The maximums and minimum images are computed as follows:

˜ t), IKm (x) =min I(x,˜ t),
IKM (x) =max I(x,(4.7)

gKM (x) =max g(x, t), gmK (x) =min g(x, t),(4.8)

T K−∆T M<t≤T K

T K−∆T M<t≤T K

T K−∆T M<t≤T K

T K−∆T M<t≤T K

where ∆T M corresponds to M times of the vibration cycle time ∆T N = 1/ fN . ∆T N is twice
the frame interval of the downsampled images. The abovementioned process is executedfor 2M frames, corresponding to the M cycle time.

32 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER(5) Image accumulation
The downsampled images IKA (x) and their high-pass ﬁltered images gAK (x) after theoﬀset processing are accumulated in a certain interval of ∆T L as follows:
KK (x) =G K (x) =Z TKT K −∆T L
Z TK
T K −∆T L

IKA (x, t)dt,(4.9)

gAK (x, t)dt,(4.10)

where ∆T L is set to L times the cycle time ∆T N .
(6) Vibration region extraction
The vibration region VK (x) is extracted by thresholding the ratio of the accumulated
images of the downsampled images KK (x) to those of their high-pass ﬁltered images
G K (x), with a threshold θ2 as follows:





 1 (G K (x) > min(θ1 , θ2 · KK (x))
VK (x) = 
,



 0
(otherwise)(4.11)

where pixel x is judged as an ambiguous pixel that should not be extracted when G K (x)is less than the threshold θ1 .
The diﬀerence of vibration source localization for HFR images at a ﬁxed frame
rate in [12] is that “(2) downsampling with frame interpolation” is added and the pixellevel bandpass ﬁlters in [12] are replaced with “(3) pixel-level high-pass ﬁlters”; the othersub-processes are similar to those in the previous method. Our algorithm can reduce its
computational complexity and required memory size in vibration source localization, corresponding to the target’s vibration frequency. It can accelerate the speed of pixel-level
digital ﬁlters by introducing, pixel-level digital ﬁlters, downsampling with frame interpolation, in which the Nyquist frame rate of downsampled images is tracked to the target’s
vibration frequency using high-pass ﬁlters, whose ﬁlter order is half that of bandpass

33

4.4 THE NUMBER OF CALCULATIONS AND COMPUTATIONAL COMPLEXITIES

Table 4.1: Computational complexities per pixel per unit time (R = f0 / fN )our proposed method
ADD / SUB
MUL
COMP
(1) image acquisition–––(2) down sampling
2 fN
4 fN
–
(3) digital ﬁlter
4P fN(4P + 2) fN–(4) oﬀset processing
4 fN–
16M fN(5) accumulation
8L fN
–
–
(6) vibration extraction
–
2 fN
4 fN
total
(8L + 4P + 6) fN (4P + 8) fN (16M + 4) fN(1) image acquisition(2) down sampling(3) digital ﬁlter(4) oﬀset processing(5) accumulation(6) vibration extraction
total

previous method
ADD / SUB
MUL
–
–
–
–
4PR fN(4P + 1)R fN
2R fN–
2LR2 fN––
R fN
2(2LR + (4P + 2)R) fN (4P + 2)R fN

COMP
–
–
–
4MR2 fN–
2R fN(4MR2 + 2R) fN

ﬁlters.
Here, we will discuss vibration source localization with pixel-level digital ﬁlters
when the translation speed of the target to be observed is not much larger than its vibration
frequency. We assume that the number of cycle times when the time-varying brightness
signal related to the target’s vibration is observable at the same pixel, and N p = A · fN /v
is suﬃciently large for computing with digital ﬁlters; A, fN , and v are the size, vibration
frequency, and translation speed of the target to be observed, respectively.

4.4

The number of calculations and computational complexities
Next, we will evaluate the number of calculations for adding and subtracting, multi-

plication, and comparison operation in this proposed method if we track vibration source
for target frequency fN (= f0 /R), which is 1/R of frame rate f0 . Table 4.4 shows the34 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

summary of the number of calculations which we need in 1 pixel per time unit. Here,
we remove the execution time regarding (1) Image acquisition. Comparing the previous
method, we also wrote the number of calculations with the previous method [12] is vibration source tracking for a ﬁxed frame rate. In the case of the proposed method, thenumber of calculations in (2) Downsampling using frame interpolation takes 1 adding
and subtracting and two multiplication every time we generate 1 frame of interpolated
image.
The number of calculations of ﬁlter order P of (3) processing high pass ﬁlter need
2P adding and subtracting and 2P + 1 multiplication. (4) Oﬀset processing with maxmin diﬀerence in M period of target frequency needs 2 adding and subtracting and 8M
comparison operation. (5) Image accumulation in L period of target frequency needs 4L
adding and subtracting. Finally, (6) Vibration source localization needs 1 multiplicationand 2 comparison operation. When we calculate the number of calculations per unit timeat 1 pixel is 2 fN times of that number. So, Total of those are number of adding and subfi
fi
tracting Dadd
= (8L + 4P + 6) fN , number of multiplication Dmul
= (4P + 8) fN , and numberfiof comparison operation Dcomp= (16M+4) fN . On the while, in the conventional vibration

source tracking method, the downsampling process of (2) is not used for comparison, but
the frame of (3) is used instead of a high-pass ﬁlter with ﬁlter order P for the input imageat frame rate f0 . A bandpass ﬁlter with ﬁlter order 2P with level performance is directly
applied, and the processing from (4) to (6) is the same as the proposed method. For eachframe of the input image, adding and subtracting is performed 4P, and multiplication is
perfomed 4P + 1 for the 2P order bandpass ﬁlter processing in (3). For the M period of
the vibration in (4), oﬀset processing using maximum/minimum value operation has 2P
adding and subtracting and 4M f0 / fN = 4MR multiplication. It is necessary to perform
adding and subtracting 2L f0 / f M = 4LR times in the integrated image calculation for L cycles of the vibration in (5), one multiplication and two comparison operations in vibration
area extraction in (6), The number of calculations per pixel in the unit time is f0 = R fN

35

Q ratio

4.4 THE NUMBER OF CALCULATIONS AND COMPUTATIONAL COMPLEXITIES

Qadd

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

100

200

300

400

500

600

700

Qmul

800

900

Qcomp

1000

fN [Hz]Figure 4.3: Relationship between computational complexities and target vibration
frequency (P = 3, M = 3, L = 6)times these, and if expressed in the form of R and fN , as shown in Table 4.4. Total numbers of calculations for the previous method are the number of adding and subtracting
pr
pr
Dadd
= (2LR2 + (4P + 2)R) fN , number of multiplication Dmul
= (4P + 2)R fN , and number
prof comparison operation Dcomp= (4MR2 + 2R) fN . The ratio of adding and subtracting

Qadd , the ratio of multiplication Qmul l and the ratio of multiplication Qcomp , per time unit
between the proposed method and the previous method are given as follows:

Qadd =

fi
Dadd
pr
Dadd

=

Dfi
4L + 2P + 3
2P + 4, Qmul = mul,pr =R((2P + 1) + LR)
Dmul R(2P + 1)
fi
Dcomp
8M + 2
Qcomp = pr =Dcomp R(2MR + 1)

(4.12)

When the target frequency fN is satisﬁed with R ≥ 2 which means target frequency
is lower than Nyquist frequency f0 /2, each of Qadd , Qmul and Qcomp are (4L+2P+3)/(4L+
4P + 2), (P + 2)/(2P + 1), 1 and they are max value when R = 2. These equations indicatethat the proposed method can be faster than the previous method if the target frequency
is less than Nyquist frequency because all of ratio are less than 1 when ﬁlter order P is

36 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

larger than 1 without considering the value of L and M.
For instance, Figure 4.3 shows the relationship between the number of calculations
ratio in Qadd , Qmul , Qcomp , and target frequency when parameters set as the ﬁlter orderP = 3, period of oﬀset processing with max-min diﬀerence M = 3, the period of image
accumulation L = 6 for input framerate f0 = 2000 fps. These parameters are following
the experiment’s setting in section 4.7; the ratio in Qadd , Qmul , and Qcomp can be small andsmall when target frequency is much smaller than input framerate. For example, when
we set target frequency fN = 200 Hz as R = 10.
The ratio of adding and subtracting Qadd , the ratio of multiplication Qmul l and theratio of multiplication Qcomp are respectively 0.049, 0.143, 0.043. As a result, we can see
the proposed algorithm is faster than the previous algorithm because we can reduce thenumber of adding and subtracting, number of multiplication and number of comparison
operation at 1 pixel per time unit.

4.5

Filter performance
When we want to speed up the function of vibration source localization, we needto realize much less order of pixel-level digital ﬁlter than one of the previous. At thesame time, we maintain ﬁlter performance as well as previous. Much research has the
purpose that it can reduce the noise for the acoustic signal and so on, while it maintains
phase characteristics and restrains the ripple of the signal. However, our proposed algorithm regarding vibration source localization has the purpose of whether it can detect
the pixel, which has a distinct frequency component for image intensity. That’s why
it is not such a severe problem to ripple in passband and stopband in a ﬁlter and nonlinear phase response with IIR ﬁlter because these image intensity’s phase periodicallychanges at diﬀerent timing pixel by pixel in general. On the other hand, setting narrowtransition bandwidth that we can separate passband with stopband is essential for ﬁlter

37

4.5 FILTER PERFORMANCE

gain

transient band
stopbandtransient band
passbandstopband

1

−

−

+

+

frequency

Figure 4.4: Transient band in digital ﬁlter.
performance. We need to localize vibration source at pixel level from the point of viewthat we can restrain misdetection for frequency component except for passband aroundtarget frequency fN . After that, we decide to consist of a pixel-level digital ﬁlter for theelliptic ﬁlter. We can design narrow transition bandwidth between passband and stopband
without compromising the ripple of passband and stopband.
As Figure 4.4 shows, we will evaluate ﬁlter performance by calculating transition
width ∆ f = ( f st+ − f st− ) − a fN after calculating stopband f < f st− , f > f st+ ( fN− = (1 −a/2) fN , fN+ = (1 + a/2) fN , a = 0.1) when we design elliptic ﬁlter which pass band fN− <f < fN+ around target frequency fN = 200 Hz and we set the parameters are maximum
passband ripple A p = 1 dB and stopband gain A st = 40 dB, These parameters are almost
same with experiment’s parameters which we set in section 4.7.
Figure 4.5 shows the gain characteristic of (a) case of using ecliptic high pass ﬁlter
whose order P = 2, 3, 5, 10 in proposed method and (b) case of using ecliptic bandpass
ﬁlter whose order is 2P in the previous method for input framerate f0 = 2000 fps. It
indicates the gain characteristic of frequency, which includes over Nyquist frequency(200 Hz) in downsampling image at 400 fps in (a) without considering double image error38 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

gain [dB]

400Hz (P=2)
2000Hz (P=2)

400Hz (P=3)
2000Hz (P=3)

400Hz (P=5)
2000Hz (P=5)

400Hz (P=10)
2000Hz (P=10)

0
-10
-20
-30
-40
-50
-60
-70
-80
-90
100

120

140

160

180

200

220

240

260

280

300

frequency [Hz]Figure 4.5: Gain properties of elliptic ﬁlters (P = 2, 3, 5, and 10).

transition bandwidth [Hz]previous (P=2)
proposed (P=2)

previous (P=3)
proposed (P=3)

previous (P=5)
proposed (P=5)

previous (P=10)
proposed (P=10)

1400

1800

140

120
100
80
60
40
20
0
400

600

800

1000

1200

1600

2000

input frame rate [fps]Figure 4.6: Relationship between input frame rate f0 and transition bandwidth ∆ f
and eﬀect because of downsampling with frame interpolation. Figure 4.6 indicates therelationship between framerate f0 = and transition width ∆ f when we change frameratefrom f0 = 420 fps to f0 = 2000 fps. In the case of ﬁlter order P = 3, transition with of (a)is constantly 27.9 Hz regardless of input framerate, that of (b) in input framerate f0 = 500,
1000, 1500, 2000 fps is 27.2, 28.0, 28.2, 28.3 Hz. So the results of (a) are almost the samewith (b).
Next Figure 4.7 shows the relationship between ﬁlter order P(= 1, · · · , 10) and transition width. It indicates transition width of (a) and (b) when we set stopband gain A st 40,
60, 80, 100 dB. If ﬁlter order P is larger and larger and input framerate f0 = 2000 fps,

39

4.5 FILTER PERFORMANCE

previous (40dB)
proposed (40dB)

transition bandwidth[Hz]400

previous (60dB)
proposed (60dB)

previous (80dB)
proposed (80dB)

previous (100dB)
proposed (100dB)

350
300
250
200

150
100
50
0
1

2

3

4

5

6

7

8

9

10

filter orderFigure 4.7: Relationship between ﬁlter order P and transition bandwidth ∆ f .
previous (P=2)
proposed (P=2)

0.7

previous (P=3)
proposed (P=3)

previous (P=5)
proposed (P=5)

previous (P=10)
proposed (P=10)

0.6

ratio

0.5
0.4
0.3
0.2
0.1
0

0

100

200

300

400

500

600

700

800

fN[Hz]

Figure 4.8: Relationship between target frequency fN and the relative ratio of transition bandwidth to target frequency ∆ f / fN
transition width of both (a) and (b) is smaller and smaller, and that of (a) is nearly corresponding with (b). For example, the transition width in A st = 40 dB are corresponding
with(a) 108.9 Hz, (b) 118.7 Hz in P = 2, (a) 27.9 Hz, (b) 28.3 Hz in P = 3, (a) 4.4 Hz,(b) 4.4 Hz in P = 5 and (a) 0.1 Hz, (b) 0.1 Hz in P = 10. If a restriction of stopband gain
is severe, the ﬁlter order needs large order to realize designing a desirable ﬁlter, which isnarrow transition width from these graphs.
Figure 4.8 indicates the relationship between target frequency fN and the relative

40 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

ratio of transition width ∆ f / fN when we set the parameter of target frequency fN from
20 Hz to 800 Hz in input framerate f0 = 2000 fps. When the ﬁlter order P is 2, therelative ratio of (a) is a little bit diﬀerent from (b) regarding target frequency. On theother hand, when the ﬁlter order P is larger than 3, the relative ratio of (a) is almost thesame with (b). For example, when the ﬁlter order P is 3, a relative ratio of (a) is 0.139
regardless of target frequency and relative ratio of (b) in the target frequency f0 = 20,
80, 200, 800 Hz is 0.145, 0.143, 0.142, 0.136. So the relative ratio of (a) is also almost
the same with (b) in this case. As a result, when we design the digital ﬁlter which is thecenter of target frequency, we can see downsampling by frame interpolation in (a) are notalmost eﬀective for the ﬁlter performance regardless of the target frequency. That’s why
we demonstrate the algorithm which is vibration source tracking with pixel-level digitalﬁlter can speed up while the ﬁlter characteristic that we need to track vibration source
maintains. This is because that the transition width especially is not so diﬀerent from thecase that we design the digital ﬁlter whose order is the same with the previous method forhigh framerate input images. It is regarding the proposed method in which we calculate
pixel-level digital ﬁlter for the downsampling image as well as the previous method.

4.6 Execution time
Next, we evaluated the execution time after we implemented the proposed algorithm
on the PC. We used a PC with an HP Z440 Workstation (Intel Xeon E5-1650v4 processer@ (3.60 GHz CPU, 6cores, 15MB, 2.4 GHz), Intel C612 chipset, 32GB RAM, and Windows 7 Enterprise 64-bit OS. Moreover, we use OpenMP [95] which makes multi-threadfor parallel processing to speed up when we implement the algorithm. We evaluated theexecution time for the proposed method, when we set the parameters as follow: vibration
frequency fN = 200 Hz, period of oﬀset processing with max-min diﬀerence M = 3 and a
period of image accumulation L = 3 as well as the parameters we set for oﬄine tracking

41

4.6 EXECUTION TIME

Table 4.2: Execution times and image sizes (unit [s]).

proposed
256
×256

previous(500 fps)
previous(1000 fps)
previous(2000 fps)
proposed

512
×
512

previous(500 fps)
previous(1000 fps)
previous(2000 fps)
proposed

1024
×
1024

previous(500 fps)
previous(1000 fps)
previous(2000 fps)
proposed

2048
×
2048

previous(500 fps)
previous(1000 fps)
previous(2000 fps)

1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr
1 thr
10 thr(1) imageacquisition(2) down
sampling(3) digitalﬁlter(4) oﬀset
processing

—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—

0.046
0.009
—
—
—
—
—
—
0.206
0.033
—
—
—
—
—
—
0.791
0.123
—
—
—
—
—
—
3.030
0.475
—
—
—
—
—
—

0.391
0.055
3.894
0.347
7.801
0.672
15.836
1.370
2.480
0.367
15.8931.476
30.3292.894
61.814
5.887
9.130
1.859
61.4845.869
122.571
11.621
175.137
23.44735.461
8.035
175.71722.867462.54745.719694.673
91.720

0.346
0.055
0.748
0.107
2.646
0.382
13.994
1.773
1.671
0.253
3.759
0.559
13.641
1.783
58.702
9.055
8.956
1.461
20.3982.664
67.995
7.863
199.587
36.78734.2785.356
58.77410.038258.16230.040791.742143.438(6) vibration(5)
regionaccumulation
extraction
0.423
0.019
0.076
0.004
1.238
0.028
0.181
0.006
5.654
0.060
0.710
0.013
25.179
0.128
3.284
0.029
2.115
0.089
0.402
0.016
6.581
0.121
0.832
0.023
24.167
0.231
3.646
0.050
179.1960.494
25.4830.105
13.800
0.338
1.910
0.067
31.271
0.456
3.660
0.089
108.2380.914
14.7870.177
524.9621.403
93.8660.370
53.336
1.314
7.476
0.270
91.4181.354
14.1590.335
411.075
3.495
56.7490.661
2067.934
5.477
361.117
1.361

total
1.226
0.199
5.909
0.640
16.162
1.777
55.136
6.456
6.561
1.072
26.3542.889
68.3688.373
300.206
40.52833.0145.419
113.60912.282299.718
34.448901.089
154.470
127.419
21.613327.263
47.4001135.279
133.1693559.826
597.636experiment for multi-copter as we will explain about it in chapter 4.7.
Table 4.2 shows the breakdown of execution time for input image as 256×256,
512×512, 1024×1024 and 2048×2048 in input framerate f0 = 2000 fps in a second’stime. Table 4.2 indicates the intermediate value of execution time when we measure theexecution time at 10 times trials with same processing. It also indicates breakdown ofexecution time regarding case of ﬁxed framerate ( f0 = 500, 1000, 2000 fps) to compare
ﬁxed framerate. The execution time in both of thread 1 and thread 10 are increasing if theimage size is increasing. In other words, we can clearly see the tendency when execution
time in previous method is longer than that in proposed method especially in case ofthat input framerate is so high. For example, if it is 512×512 image at input framerate
f0 = 500, 1000, 2000 fps, execution time in previous method is longer than proposed

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

2048x2048

512x5121024x1024

256x2562048x2048

512x5121024x1024

256x256(6) vibration region extraction
2048x2048

512x5121024x1024

256x2562048x2048

512x512(5) accumulation
1024x1024

256x2562048x2048

512x5121024x1024

256x256(4) offset processing
2048x2048

512x5121024x1024

256x2562048x2048(3) digital filter
512x512256x2562048x2048

80%

512x51290%

256x256100%

1024x1024(2) down sampling

1024x1024

42 CHAPTER 4.

70%
60%
50%
40%

30%
20%
10%
0%

10 thread
1 thread
proposed

10 thread
1 threadprevious (500)

10 thread
1 threadprevious (1k)

10 thread
1 threadprevious (2k)

Figure 4.9: Breakdown percentages of execution times.
method as follow : 4.0, 10.4, 45.8 times in thread 1 and 2.7, 7.8, 37.8 times in thread 10.
Both method’s execution time of proposed and previous in thread 10 are shorter than thatin thread 1 regardless of image size as its ratio is about execution time in previous method
is longer than proposed method as follow : 4.0, 10.4, 45.8 times in thread 1 and 2.7,
7.8, 37.8 times in thread 10. Both method’s execution time of proposed and previous inthread 10 are shorter than that in thread 1 regardless of image size as its ratio is about 16 – 16
by using multi thread processing. Figure 4.9 indicates occupancy of each processing forexecution time. And both method’s occupancy is mainly consisted of (3) digital ﬁlter, (4)oﬀset processing and (5) Image accumulation in more than 90% regardless of image sizes
and number of thread. In this case, we can see (2) down sampling which we introduced
in proposed method is not interfere with other part of processing.
Next, Figure A.2 shows the relationship between execution time for input image
in a second’s time in both method and image size when we make number of pixel Mchange from 20 to 1000 by every aproximately10 pixels in case of image size is M × M
pixel. The parameters show as follow: input framerate f0 = 2000 fps, number of thread
10 and ﬁlter order P =3. Execution time for input image in a second’s time in Figure
A.2 is within 1 second. In other words, maximum size of image which we can execute

43

ratio

4.6 EXECUTION TIME

proposed

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

1

2

3

4

5

previous (500)

6

7

8

previous (1k)910previous (2k)1112threadFigure 4.10: Relationship between execution times and image sizes.
on real time using thread 10 is 550×550 pixels in proposed method while that in previous
method at input framerate f0 = 500,1000,2000 fps is 340×340,240×240,130×130 pixels
respectively. So the number of pixels in previous method we can process on real time is0.40, 0.20, 0.06 times as much as that in proposed method. As you can see these results, it
is possible for proposed algorithm to realize vibration source localization for much higher
space resolution than previous algorithm.
Figure 4.11 indicates relationship between execution time and number of threadwhen we set number of threads as 1–12 in proposed algorithm (P=3) for 512×512 input image in a second’s time at input framerate f0 = 2000 fps. Here, we show relativechange of execution time in previous method and proposed method for number of threads
as the execution time in thread 1 is standard. The execution times in thread 1 are 6.561
s in proposed algorithm and 23.314, 60.761, 268.990 s in previous algorithm ( f0 = 500,
1000, 2000 fps). When number of threads is larger than 10, execution regarding numberof threads is almost same in both algorithm while execution time in previous and proposed algorithm is decreasing when number of threads is increasing. The execution times
in thread 10 are 0.971 s in proposed algorithm and 2.458, 7.232, 37.555 s in previous
algorithm ( f0 = 500, 1000, 2000 fps). Comparing with thread 1, execution time in previous algorithm is 0.15 times and execution times in previous algorithm( f0 = 500, 1000,

44 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

proposed3

previous (500)previous (1k)previous (2k)

2.5

ratio

2
1.5
1
0.5

0
1

2

3

4

5

6

7

8

9

10

filter order P

Figure 4.11: Relationship between execution times and number of threads.

processing time[s]proposal

previous (500)previous (1k)previous (2k)

10
9
8
7
6
5
4
3
2
1
0

0

100

200

300

400

500

600

700

800

900

1000

input size M [pixel]

Figure 4.12: Relationship between execution times and input frame rate.
2000 fps) are respectively 0.12, 0.14, 0.15 times, so the execution time in thread 10 speed
up around 7 times faster than thread 1. These execution times are depending on numberof PC processor cores which we evaluated, and it is corresponding to implementation ofvibration source tracking algorithm with multi thread by synchronizing Windows OS withother standard software for processor Intel Xeon E5-1650v4 which provide 12 processing
thread with 6 physical cores.
Figure 4.12 indicates execution time in proposed and previous method (P =3) for
512×512 input image in a second’s time in thread 10 when input framerate f0 changes
between 400 and 2000 fps. The execution time in proposed algorithm is always 0.971 s

45

4.6 EXECUTION TIME

proposed45

previous

execution time [s]

40
35
30

25
20

15
10
5
0
400

600

800

1000

1200

1400

1600

1800

2000

input frame rate [fps]Figure 4.13: Relationship between execution times and ﬁlter orders.
by introducing down sampling processing with frame interpolation regardless of input
framerate. On the other hand, the execution time in previous algorithm ( f0 = 500, 1000,
2000 fps) is increasing with increasing input framerate as we explained before and indicates 2.707, 7.866, 41.150 s respectively. For example, prosed method is over 40 times
faster than previous one in case of f0 = 2000 fps. The proposed method dramatically can
expect to speed up for previous method when input framerate f0 is much higher than target frequency fN = 200Hz because computational complexity is decided by dependenceof the target frequency fN = 200Hz.
Last, Figure 4.13 indicates the relationship between execution time and ﬁlter order
P in both previous and proposed method for 512×512 input image in a second’s time at
input framerate f0 = 2000 fps, when we set ﬁlter order P = 1–10. In Figure 4.13, we
can show relative change of execution time in previous method and proposed method forﬁlter order as the execution time in ﬁlter order P =3 is standard. The execution times inproposed algorithm and in previous algorithm ( f0 = 500, 1000, 2000 fps) are increasing
as ﬁlter order P is increasing. In addition, execution times in ﬁlter order P =10 is 2.7
times in proposed method and 2.3, 1.8, 1.3 times longer than that in previous algorithm
in P =3. Moreover, the eﬀect of execution time in previous algorithm for ﬁlter orderchanging tends to be small when input framerate f0 is higher than target frequency fN .

46 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

Figure 4.14: Experimental environment when video-shooting a ﬁxed multi-copter.
This is because that computational complexity of (4) oﬀset processing and (5) image
accumulation which need to process the number of periods for vibration frequency as
number of frame in a vibration period is increasing with increasing input framerate isincreasing regarding that ﬁlter order P for (2) digital ﬁlter in previous method can beﬁxed value.

4.7

Experiments for a multi-copter

4.7.1

Video shooting a ﬁxed multi-copter with various exposure timeFirstly, we checked wheter we could extract region in which multi-copter’s pro-peller rotate as vibrated region in proposed algorithm as well as previous vibration source
localization. Considering about that there are double image errors which occurs in casethat exposure time is longer than the period of framerate in linear interpolated image afterrecording 512×512 image of 8 bit color for a ﬁxed multi-copter and fans at 2000 fps, when
we set exposure time as 0.5, 0.2, 0.1 m. We selected to use multi-copter (RC EYE One
Extreme by CEI Conrad company) which is size of object is 225×225×80 mm and have
4 two-blade propellers whose each size is 138 mm for evaluation. Figure 4.14 indicates47

4.7 EXPERIMENTS FOR A MULTI-COPTER

input images

t=0.000s

t=0.250s

frame interpolated images

t=0.500s

t=0.000s

t=0.250st=0.500s(a) exposure time 0.5 ms
input images

t=0.000s

t=0.250s

frame interpolated images

t=0.500s

t=0.000s

t=0.250st=0.500s(b) exposure time 0.2 ms
frame interpolated images

input images

t=0.000s

t=0.250st=0.500s

t=0.000s

t=0.250st=0.500s(c) exposure time 0.1 m
Figure 4.15: Input images when video-shooting a ﬁxed multi-copter.
these experimental environment. We set the distance as 4.5 m from camera to multicopter and a pixel is corresponding to 1.0 mm for multi-copter’s position and 512×512
image is corresponding to 512×512 mm. In addition, we ﬁxed the body of multi-copterat an angle in order to monitor all of motion of 4 propellers. And we also ﬁxed 2 twoblade fans whose rotation speed is set as 250, 60 rotation per second respectively next tomulti-copter as we also can see background of trees on sun light. We make multi-copter’s
propeller rotating at 100 rps in all setting in this section. We set parameters in proposed
algorithm as follow: vibration frequency fN = 200 Hz, ﬁlter order P=3, period of oﬀset
processing with max-min diﬀerence M =3, period of image accumulation L =6.
Figure 4.15 indicates input frames and frame interpolated images’ group in proposed method when we set exposure time (a) 0.5 ms, (b) 0.2 ms, (c) 0.1 ms, frame inter-48 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

proposal

previous(500)previous(1k)previous(2k)(a) exposure time 0.5 ms

proposal

previous(500)previous(1k)previous(2k)(b) exposure time 0.2 ms

proposal

previous(500)previous(1k)previous(2k)(c) exposure time 0.1 ms
Figure 4.16: Extracted vibration features when video-shooting a ﬁxed multi-copter.

Figure 4.17: Experimental environment for ﬂying multi-copter.
polated images can be double image error when exposure time is longer especially while
motion blur in propeller’s region can be reduced and input image turns dark all over theimage.
Figure 4.16 shows vibration region by using proposed method as we compare vibration region by implementing previous algorithm ( f0 = 500, 1000, 2000 fps) in the case4.7 EXPERIMENTS FOR A MULTI-COPTER

49

of exposure time (a) 0.5 ms, (b) 0.2 ms, (c) 0.1 ms. We can detect areas in which there
are multi-copter’s propeller rotating at 100 rotation per second as vibration region while
we don’t detect areas in which there are ﬁxed 2 two-blade fans whose rotation speed is
set as 250, 60 rps as vibration region in both methods. So this tendency is not depending
on exposure time. Furthermore, in terms of the number of detected pixels, when input
framerate f0 = 2000 is in the ﬁrst case by 100 and 500 fps for previous method in Figure 4.16. This is because that there are some pixels whose diﬀerence are small between
max and min, and number of sampling is not enough to extract actual intensity change incase of low framerate. So the number of pixels detected in equation (4.11) by calculating
30 frames, 15 frames and 8 frames ( f0 = 500, 1000, 2000 fps) regarding maximum andminimum brightness for 15 ms ( fN = 200 Hz, M = 3), which is corresponding to equation (4.7). On the other hand, number of detected pixels in proposed algorithm in equation(4.11) is decreasing because Nyquist frequency is corresponding to target frequency andnumber of pixel which we can detect maximum and minimum brightness as its phase isalmost inversely at 180 degrees is increasing and number of pixel whose diﬀerence between max and min is much larger than previous method regardless of input framerate
f0 = 400 fps. As a result of these considerations, we can realize vibration source tracking
in proposed algorithm performs as well as that in previous algorithm regardless of double
image errors when the exposure time is small.

4.7.2

Video shooting a ﬂying multi-copter
Next, we will show results of vibration source tracking for 512×512 image of 8 bitcolor for a ﬂying multi-copter in outdoor at 2000 fps as background is walking humans
and trees and so on when we apply the proposed algorithm. We used multi-copter whichis same one in chapter 4.7.1. Figure 4.21 indicates the experimental environment. Here,
we can monitor brightness change in enough period at the same pixel in this experimentbecause of parameter N p = 35.4 as followed two-blade propeller’s size A is 138 mm,

50 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

target frequency fN = 200 Hz, maximum speed v = 0.78 m/s. We recorded the experimental images at 2 patterns as (1) case with large magniﬁcation and (2) case with wide
angle. We set the distance as 24.0 m from camera to ﬂying multi-copter and 512×512
image is corresponding to 3072×3072 mm in (1) case and 7680×7680 mm in (2) case formulti-copter’s position and a pixel is corresponding to 6.0 mm in (1) case and 15.0 mm
in (2) case. In addition, we used the parameters are same with previous indoor experiment. Figure 4.18 indicates time sequential images of (a) input images with tracked
positions, (b) accumulated images of input KK (x), (c) pixel-wise ﬁltered images G K (x),(d) vibration features extracted by our proposed algorithm VK (x), (e) vibration regions extracted by previous algorithm ( f0 = 2000 fps) and (f) ROI input images for tracked regions
whose gravity position is center of this image, when video-shooting a ﬂying drone withlarge magniﬁcation. The rectangular region in (a) is corresponding to magniﬁed 50×50
ROI image of (f). We can extract the vibration region which is corresponding to regionof multi-copter’s propeller rotating in proposed algorithm as well as previous algorithm
even though multi-copter moved around areas in full image. In addition, we can specify
the vibration source region which is same as position of ﬂying multi-copter because there
is always multi-copter in center of ROI image of (f).
Figure 4.22 indicates time sequential image of (a) input images with tracked positions, (b) accumulated images of input, (c) pixel-wise ﬁltered images, (d) vibration features extracted by our proposed algorithm and (e) previous algorithm ( f0 = 2000 fps) and(e) magniﬁed images for tracked regions when video-shooting a ﬂying drone with a wider
view angle. Even though multi-coper’s image is very low resolution’s image which corresponds to 30×30 image as Figure 4.22 (f) shows it, comparing with zoom image whichcorresponds to about 50×50 pixels’region, we can realize pixel level vibration region extraction for region of multi-copter’s propeller in proposed algorithm as well as that inprevious algorithm regardless of eﬀect for noise which is moving humans and somethingaround natural tree in background. We can track multi-copter as vibration region because

4.7 EXPERIMENTS FOR A MULTI-COPTER

51

there is multi-copter in center of ROI image in almost all time although there are somecases that the position of multi-copter is a little bit diﬀerent from center of the image due
to some occlusion which the propeller sometimes overlaps multi-copter itself.
Figure 4.23 shows (a) time sequential data of gravity position and number of pixels regarding vibration region for 10 seconds and (b) xy trajectory of extracted vibration
region when a ﬂying drone was video-shot with a wider view angle. In order to comparewith the results of extracted region in previous method, it also is shown in Figure 4.23
Moreover trajectory in Figure 4.23 (b) is plotted on input image in time t =0 in Figure
4.22 (a). The gravity position of vibration region in proposed method nearly corresponds
to that in previous method as it is same with up, down, left and right of movement ofmulti-copter even though number of pixel in extracted vibration region change from 0 to
54 signiﬁcantly when multi-copter moved in front of walking human and trees as background. Furthermore, we also can realize robust vibration source tracking with pixel level
digital ﬁlter for ﬂying multi-copter though there are some occlusions which the propeller
sometimes overlaps multi-copter itself and some cases that the gravity position of multicopter is a little bit diﬀerent from center of ROI image.

4.7.3

Real-time multi-copter detection
We show the experimental results of real-time vibration source detection when amulti-copter moved against a background with trees and moving persons. The multicopter ﬂew to the right and left directions above many people with quick arm and leg
movements. The 512×512 input images were captured at 2000 fps with 0.5 ms exposure.
We used multi-copter which is same one in chapter 4.7.1. The ﬂapping frequency of each
propeller varied by approximately 95 rps. Figure 4.21 shows the experimental environment. The multi-copter ﬂew at a distance of approximately 17.4 m from the camera. The
measurement area was 5.24×5.24 m for 512×512 pixels at a distance of 17.4 m from thecamera, where one pixel corresponds to 10.2 mm. The target vibration frequency was set

52 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

to fN = 190 Hz and the pass bandwidth in the digital ﬁlters was 20 Hz so that the pixellevel digital ﬁlters can pass the frequency components of time-varying brightness signals
in the range of 180–200 Hz. The parameters were set to M = 3, L = 6, and P = 3. Ouralgorithm was executed in real time with multithread processing of 10 threads.
Figure 4.22 (a)–(c) shows the input images, the accumulated images of input andhigh-pass ﬁltered images. The images were taken at intervals of 1 s for t = 0–5 s. Figure 4.22 (d)–(e) shows the vibration pixels extracted by our algorithm and the previous
algorithm in [12] when f0 = 2000 fps. Figure 4.22 (f) shows the magniﬁed images of
80×80 pixels around the vibration pixels extracted by our algorithm, corresponding to
the red-line windows illustrated in (a). When the multi-copter ﬂew in the right and left
directions against the cluttered background with trees and moving persons, our algorithm
extracted certain pixels around the propellers of the multi-copter in real time. We canconﬁrm that the regions of interest illustrated in Figure 4.22 (f) included the multi-copter
at all the times, whereas their space resolution was too low for precise appearance-based
recognition.
Figure 4.23 illustrates (a) the graphs that show changes in the x- and y-coordinate
values of the averaged positions of the extracted pixels and the number of extracted pixels
for 5 s; and (b) the xy trajectory for 5 s was plotted over the input image at t = 0.
Corresponding to the left and right movements of the ﬂying multi-copter, the xy trajectoryof the averaged positions of the extracted pixels were robust without any disturbance fromthe cluttered background of the outdoor location.

4.8 Concluding remarks
In this chapter, we demonstrated vibration source tracking for object’s frequency is
200 Hz with input images at 2000 fps in proposed algorithm is over 10 times faster thanthat in previous method after we propose a fast algorithm which reduce computational

4.8 CONCLUDING REMARKS

53

complexity per unit time by controlling sampling period for digital ﬁlter with downsampled HFR images whose frame rates were virtually adjusted to twice the target vibration
frequency with frame interpolation, while locking in vibration frequency for tracking object, we check the number of calculations and ﬁlter performance and we evaluate execution time on PC with CPU regarding vibration source tracking for high framerate video
images with pixel level digital ﬁlter. We conﬁrmed that the execution time of our algorithm was 1/10 or less than that of the previous algorithm when localizing a vibration
source using HFR images captured at 2000 fps. The eﬀectiveness of this porposed algorithm was veriﬁed by the experimental results for HFR oﬄine video analysis and real-time
detection for multi-copter, in which a multi-copter with propellers rotating at 95 rps thatﬂew in a cluttered outdoor location was simultaneously localized using 512×512 images
at 2000 fps. Furthermore, in this section, we focused on speeding up for pixel level vibration source localization which can extract frequency information in image pixel by pixel
and we discussed with the case that frequency of multi-copter’s propeller is obviously
higher than frequency of background and the other object’s motion. However, robust
recognition method which integrated dynamics characteristic for acoustic frequency levelwith space appearance characteristic can be meaningful if we consider about the case that
there are various kinds of dynamic background and multiple similar objects in a pointof view in actual scene. Moreover, vibration source localization for downsampled image
which controlled locking in the sampling period for peak frequency can also be valid bydetecting peak vibration frequency based on intensity change for acoustic signal and image because the number of rotations of propeller and so on is not constant for conditionof multi-copter ﬂying and depending on various kinds of that. Vibration source tracking
algorithm we proposed in this section will be able to process much more pixels than now
and we will improve this system which can track multi objects at the same time by utilizing that we can set variable framerate corresponding to various vibration frequency with
software after implementing the function of detecting vibration frequency automatically54 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

in future work. Finally, we will improve our vision-based vibration source localization
system with real-time feedback control of its frame rate so that it can robustly localize andtrack a ﬂying multi-copter when the rotation frequencies of its propellers are unknown orvary over time..

4.8 CONCLUDING REMARKS

55

Figure 4.18: (a) Input images with tracked positions, accumulated images of (b)input and (c) pixel-wise ﬁltered images, vibration features extracted by (d) our proposed algorithm and (e) previous algorithm ( f0 = 2000 fps), and (e) magniﬁed images
for tracked regions when video-shooting a ﬂying drone with large magniﬁcation.

56 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

Figure 4.19: (a) Input images with tracked positions, accumulated images of (b)input and (c) pixel-wise ﬁltered images, vibration features extracted by (d) our proposed algorithm and (e) previous algorithm ( f0 = 2000 fps), and (e) magniﬁed images
for tracked regions when video-shooting a ﬂying drone with a wider view angle.

XY

ON pixels

x- and y-cordinate [pixel]

500
400
300
200
100
0
0

1

2

3

4

5

time [s]

6

7

8

9

100
90
80
70
60
50
40
30
20
10
0

ON-pixels

57

4.8 CONCLUDING REMARKS

10(a) x- and y-coordinates and number of pixels(b) xy trajectory
Figure 4.20: xy trajectory of extracted vibration region when a ﬂying drone wasvideo-shot with a wider view angle.

58 CHAPTER 4.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

Figure 4.21: Experimental environment for ﬂying multi-copter.

4.8 CONCLUDING REMARKS

59

Figure 4.22: (a) Input images; accumulated images of (b) input; (c) high-pass ﬁlteredimages; vibration pixels extracted by (d) our algorithm; and (e) previous algorithm( f0 = 2000 fps); (f) magniﬁed images.

A FAST VIBRATION SOURCE TRACKING ALGORITHM USING PIXEL-LEVEL DIGITAL FILTERS FOR A FLYING MULTI-COPTER

500

X

Y

200

ON pixels

400

150

300
100
200

ON-pixels

x- and y-cordinate [pixel]

60 CHAPTER 4.

50

100
0

0
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

time [s](a) x- and y-coordinates and number of pixels(b) xy trajectory
Figure 4.23: The xy trajectory of extracted vibration region when the video of theﬂying drone was shot with a wider angle.

Chapter 5
HFR-video-based honeybee activity sensing
5.1

Introduction
Presently, the honeybee industry is of huge beneﬁt to agricultural businesses bypollinating the ﬂowers of crops, vegetables, and fruits as well as by producing and storinghoney; the honeybee is recognized as one of the most eﬃcient insect pollinators that
can transport and store pollen [96]. Colony collapse disorder (CCD) is the phenomenon
where the majority of worker bees in a colony disappear rapidly, leaving behind the queen.
The number of disappearances of beehives has increased drastically in the 2000s [97];
CCD has become a serious social problem that causes signiﬁcant economic losses toagricultural business because honeybees pollinate many agricultural crops. Whereas thecauses of CCD are not perfectly clariﬁed yet, various possible causes such as plagues
and viruses, malnutrition, pesticides and miticides, beekeeping practices, electromagnetic
radiation, and genetically modiﬁed crops have been reported [98–100]. To obtain usefulknowledge for eﬃcient beehive management without CCD outbreak, it is necessary tomonitor colony activity to elucidate the reason for this syndrome [101].
Honeybee colonies monitoring had been attempted already more than 100 years
ago; Gates [102] have reported the temperature data of a beehive collected manually every hour over several days. The rapid progress of sensor technology in recent years has
resulted in the use of various types of electronic sensors such as (1) weight, (2) tempera61

62

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

ture, humidity, and gas, (3) sound and vibration, and (4) forager traﬃc [103] to conduct
the continuous monitoring of honeybee colonies.
Beehive weight is related with colony growth and daily colony food demand as wellas colony component such as adults, brood, and food reserves. To provide information onvarious issues such as weather eﬀects [104], colony growth and consumption [105], andhive abandonment [106], continuous weighing with a precise electronic balance has been
conducted. Temperature, gas concentration, and humidity are related with the metabolic
processes of the colony level; calculations from temperature and O2 consumption can beused to estimate the energy produced by a colony. Monitoring with these modalities has
been conducted to observe colony temperature control against decreasing ambient temperature [107], low O2 level due to reduced metabolic rate when water / energy conservation,and longevity [108], daily cycle in O2 linked to forager activity [109], and time-varying
humidity near the brood in a beehive [110].
Vibration and sound observed in a hive or its surface involve biological information that bees generate wing-ﬂapping vibrations at many frequencies in the range of 10
to 1000 Hz to communicate such as the waggle dance [111–114]. The wing-ﬂapping
vibration of a single bee was measured with laser vibrometry and analyzed with STFT
to examine it in the time and frequency domains [115]. By embedding an accelerometer in the wall of hives, Bencsik et al. [116] analyzed its vibration data with principal
components analysis to examine its vibration frequencies to the noise in swarming eventsrobustly. By conducting 270-hour monitoring, several swarming events were observedwith increased sound intensity and decreased temperature and humidity in the hives [117].
One-year sound intensity data with several frequency bands was analyzed to quantify how
the sound in healthy hives diﬀers from that in chalkbrood infected hives [118].
Forager activity [119, 120], indicating how many older bees were gathering food, is
important to monitor honeybee colony health as well as food availability / demand, andcolony age structure. Hive entrance counters such as an electrical counter with a balance5.1 INTRODUCTION63

arm [121], photoelectric sensors [122, 123], and infrared light sensors [124] have beendeveloped. Danka et al. [125] measured forager traﬃcs during blueberry pollination, andfaced a miscounting problem with photoelectric sensors because the bees tended to clusteraround the sensors; they improved the measurement accuracy by validating forager ﬂights
with a visual count. Several attempts at forager traﬃc monitoring have been reported byattaching RFID chips to bees [126, 127]; they can collect crucial information such as life
expectancy time and foraging time because bees tagged with chips are identiﬁed individually. The chips are small and weigh a few milligrams, however, it was too laborious
and expensive to attach them to numerous bees, and the hive entrance should be modiﬁed to set up RFID receivers. Several lidar-based systems [128–130] have been reported
for the large-area detection of honeybees. By detecting their wing-ﬂapping in the audio
frequency range using mechanical laser scanning, they can measure the 3-D positions of
ﬂying honeybees at a distance of hundreds of meters. However, they cannot keep track ofa speciﬁc ﬂying honeybee to monitor its ﬂying activity for a certain interval.
Activity analysis of honeybees inside beehives using automatic bee tracking and behavior labeling has been reported in [131, 132]. Many systems and studies have analyzed
the ﬂight activities of honeybees outside beehives with capacitance-sensor-based monitoring [133], wingbeat frequency analysis [134], radio-frequency identiﬁcation (RFID)
chips for honeybee counting [135], and computer-vision-based monitroing such as honeybee detection with background subtraction [136, 137], automated bee-counting system [138, 139], and three-dimensional bee tracking [140]; they used standard videos
recorded at dozens of fps, in which the wing-ﬂapping of bees were unobservable, andbecause appearances of small honeybees in images became unclear without large magniﬁcation, analyzable scenes were limited to small areas around the hive entrance. Mostof the researchers have handled honeybee monitoring in an exceedingly limited narrow
area around a beehive. Especially in vision-based analysis, there is a tradeoﬀ between themeasurement area and detection accuracy, because the size of a honeybee is a dozen of

64

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

millimeters; honeybees cannot be distinguished manually by their appearances in images
when the measurement area such as 1×1 m or more is much larger than their size.
As we introduced appearance-based tracking methods in subsection 2.2, the limitation of these methods in tracking honeybees ﬂying in a natural scene is the appearancesof honeybees change and become ambiguous in the low resolution images. By detecting
temporal periodic changes in brightness signals, HFR video analysis has been applied tovibration distribution measurement such as structural vibration analysis [141–143] andvocal fold vibration [144–146]. If honeybees with wings ﬂapping at hundreds of Hertz
could be observed in an HFR video, activity sensing of honeybees ﬂying within an areathat was a few meters in diameter around a beehive, which was much larger than the sizeof a honeybee, could be realized without attaching sensors by using pixel-level digital
signal processing; this will help to quantify the colony health before and after CCD orother serious incidents. Pixel-level digital signal processing was used for honeybee localization in [147] to extract pixels around honeybees. When a honeybee is not hovering inthe same place, however, it is diﬃcult to precisely obtain its wing-ﬂapping frequency as
its ﬂying activity. This is because one or no wave is observed in the image intensities at
pixels around fast-ﬂying honeybees.
Therefore, we propose a concept of honeybee activity sensing that can precisely
obtain the wing-ﬂapping frequencies of fast-ﬂying honeybees as well as their positions
by executing two-step short-time Fourier transforms (STFTs) for an HFR video, which
is an extension of the honeybee localization algorithm in [147]. Pixel-level STFTs are
executed at all the pixels of an HFR video as the 1st-step process for honeybee localization
in a similar manner as in [147], and pixel-level STFTs of the image intensities aroundhoneybees are re-executed at all the pixels of the tracked region-of-interest (ROI) images
as the extended 2nd-step process to precisely obtain their wing-ﬂapping frequencies. Thisis because honeybees are virtually located at the centers in the tracked ROI images, andmultiple waves can be observed for precise frequency response computation.

5.2 CONCEPT

5.2

65

Concept
Western honeybees are often used in beekeeping all over the world; the sizes ofadult worker bees are 12–14 mm, and their wing-ﬂapping frequencies when ﬂying are
170–270 Hz in the audio frequency range [148]. It is diﬃcult for a video-rate camera
operating at dozens of fps to detect such a small ﬂying honeybee under natural complex
background conditions. Regarding an HFR camera as a collection of photo sensors, thisstudy proposes a quantiﬁcation method for the ﬂying activities of multiple honeybees inimages with pixel-level digital ﬁlters. Our method is an extension of the vision-based
vibration source localization approach [12], in which a pixel-level digital ﬁlter had beenused for localization whereas no pixel-level frequency response analysis had been considered. Figure 5.1 shows the concept of HFR-video-based activity sensing that this study
presents.
Fast-ﬂying honeybees in an HFR video are localized with pixel-level digital signal processing when their appearances are poor in the low resolution and blurred images [147]. There are still many problems to precisely obtain their wing-ﬂapping frequencies as their ﬂying activities. Here, the maximum ﬂying speed of a western honeybee wasreported as 25 km/h in [149], and because this movement is so large, one or no wave isobserved at the same pixel. When a 12-mm-size honeybee is ﬂying at 25 km/h, the dwell
time in a single pixel, where its wings are detectable, is 1.7 ms in the worst case; it issmaller than one cycle time of its wing-ﬂapping. Figure 5.2 illustrates the problems inpixel-level digital signal processing for fast-ﬂying honeybees. The temporal sequencesof image intensities at the pixels to be analyzed are largely dependent on backgroundscenes. Excluding a short time for one-or-no-wave observation of their ﬂapping wings, it
is diﬃcult to use pixel-level digital signal processing to obtain precise TFRs in image intensities at pixels around fast-ﬂying honeybees. A distinct peak in the frequency response
might not be detected for one-or-no-wave observation, as illustrated in Figure 5.2(a). In
addition, the latency eﬀect in pixel-level digital signal processing [12, 13] enlarges when

66

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

a fast-ﬂying honeybee is being observed; several pixels are incorrectly detected as ghost
vibration regions even when the honeybee is not currently at these pixels owing to thelatency eﬀect in the pixel-level STFTs, as illustrated in Figure 5.2(b). A 12-mm-size honeybee at a ﬂying speed of 25 km/h moves 0.44 m during the processing time of 64 ms inexecuting STFTs with 32 samples for a 500-fps video; it is 37 times the length of its body.
Considering these issues in pixel-level digital processing, our HFR-video-based activity sensing in this study, is designed with two-step pixel-level STFTs for the precise
estimation of the time-variant ﬂapping frequency responses of individual ﬂying honeybees as well as their locations and velocities, when numerous honeybees are observed in
the ﬁeld of view of an HFR camera by executing the following processes as illustrated in
Figure 5.1.

5.3

Proposed algorithm
On the basis of the concept in Chapter 5.2, we implement a pixel-level signal pro-cessing algorithm in this Chapter, with a two-step pixel-level STFT method for HFRvideo-based honeybee activity sensing, in which the 1st pixel-level STFTs are used forvibration source localization and the 2nd pixel-level STFTs are used for precise TFR
computation. The details are provided in the following subsections.

5.3.1

Pixel-level vibration source localization
The input images of M × N pixels are acquired at time kτ (frame number k andframe cycle time τ) as I(x, y, kτ), wherein the frame rate is f0 = 1/τ. The input images
are converted to the TFR image by executing the STFT with K1 (= 2K1′ ) samples:F(x, y, t) = (F0 (x, y, t), · · · , F K1′−1 (x, y, t)),= STFT(I(x,y, t−K1′ τ), · · · , I(x,y, t+(K1′ −1)τ),

(5.1)

67

5.3 PROPOSED ALGORITHM

where K1 is the number of samples in the STFT for the input images, and F K (x, y, t)
indicates the frequency component image at a frequency of fK = f0 K/K1 (K = 0,· · ·, K1′−1).
The absolute value image of the cut-oﬀ TFR image F′k1 (x, y, t), in which frequency
components less than fk1 are cut oﬀ from F(x, y, t), is normalized by that of the DC-cutoﬀ TFR image F′k0 (x, y, t), whose cut-oﬀ frequency is fk0 as the following high-pass-ﬁlter
image:

H(x, y, t) =|F′k1 (x, y, t)||F′k0 (x, y, t)|

,

(5.2)

where ki (i = 0, 1) is an integer to set the cut-oﬀ frequency fki = f0 ki /K1 of the cut-oﬀ TFR
image as follows:
F′ki (x, y, t) = (0, · · · , 0, Fki (x, y, t), · · · , F K1′ −1 (x, y, t)).

(5.3)

The vibration candidate region is detected pixel-wise by checking the DC-cut-oﬀ
TFR image F′k0 (x, y, t) and the high-pass-ﬁlter image H(x, y, t) as follows:



′


 1 (H(x, y, t) ≥ θH , |Fk0 (x, y, t)| ≥ θF )V(x, y, t) = 
,



0
(otherwise)(5.4)

where θF is a threshold to determine whether the amplitude of the DC-cut-oﬀ TFR image
is detectable. θH is a threshold to determine whether the high-pass-ﬁlter image, whichindicates the ratio of the frequency components greater than or equal to the cut-oﬀ frequency fk1 , is suﬃciently large.

68

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

5.3.2 Individual honeybee labeling
The zeroth- and ﬁrst-moment features are calculated for M ′ N ′ subregions Γab ofm × n pixels in V(x, y, t) as
M pq (Γab , t) =a(m+1)−1
X b(n+1)−1X
x=am

x p yq · V(x, y, t) ,

(5.5)

y=bn

where M = mM ′ and N = nN ′ ; p and q are non-negative integers to determine the orderin moment feature calculation, satisfying p + q ≤ 1. Γab (a = 0,· · ·, M ′ −1, b = 0,· · ·, N ′ −1)
is expressed asΓab = {(x, y)| (am+ s, bn+t), 0 ≤ s < m, 0 ≤ t < n}.

(5.6)

To localize multiple honeybee regions in V(x, y, t), the connected components labeling process is accelerated to obtain the labeled honeybee regions with the cell-basedlabeling algorithm using the moment features for M ′ ×N ′ subregions [?]. For the labeledhoneybee regions l O(t) (l = 0,· · ·, L − 1) in V(x, y, t), the label-domain moment features
M pq (l O(t)) (p+q ≤ 1) are sequentially accumulated, concurrently with the scanning of aﬂag map Pab (t) of M ′ ×N ′ subregions:
M pq (l O(t)) =Xx p yq · V(x, y, t),(5.7)(x,y)∈l O(t)where Pab (t) (a = 0,· · ·, M ′−1, b = 0,· · ·, N ′−1) is deﬁned for each subregion Γab by checking
M00 (Γab , t) as follows: 1 (M00 (Γab , t) ≥ θP )
Pab (t) = 
,



 0
(otherwise)

(5.8)

69

5.3 PROPOSED ALGORITHM

and θP is a threshold to determine whether the subregion Γab is used for cell-based labeling.
The center positions of the labeled honeybee regions l O(t) are calculated as follows:!
M10 (l O(t)) M01 (l O(t))x(t) = ( x(t), y(t)) =,.
M00 (l O(t)) M00 (l O(t))

l

l

l

(5.9)

The detailed processes used to calculate the label-domain moment features in animage are given in [19, 20].

5.3.3

Frequency response computation for tracked honeybees
The labeled honeybee region l O(t) is deﬁned as a sub-region of R×R pixels so thatthe honeybee locates in its center, and its ROI image lI(x′ , y′ , t) of R×R pixels is selectedfrom the input image as
I(x′ , y′ , t) = I(x′ + lx(t)−R/2, y′ + ly(t)−R/2, t),l(5.10)

where (x′ , y′ ) indicates the local coordinate system in the ROI image of R×R pixels.
For the ROI images, the image intensities at pixel (x′ , y′ ) are converted to the TFR
ROI images by executing the STFT with K2 (= 2K2′ ) samples as follows:
F(x′ , y′ , t) = (lF0 (x′ , y′ , t),· · ·, lF K2′−1 (x′ , y′ , t)),l

= STFT(lI ′ (x′,y′, t−K2′ τ),· · ·, lI ′ (x′,y′, t+(K2′ −1)τ),(5.11)

where K2 is the number of samples in the STFT for the ROI images, and where l F K (x′ , y′ , t)(l = 0,· · ·, L−1) indicates the frequency component of the TFR ROI image at a frequency
of fK′ = f0 K/K2 (K = 0,· · ·, K2′ −1).
The vibration pixels in the ROI images, l V(x′ , y′ , t) (l = 0,· · ·, L − 1), are extracted
when the peak frequencies are greater than or equal to the cut-oﬀ frequency fk2 , and70

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

their TFR images at the peak frequencies are greater than or equal to a threshold θ′F todetermine whether the frequency spectrum peak is distinct:



l
′ ′
l
′ ′
′


 1 ( f pk (x , y , t) ≥ fk2 , Flk pk (x , y , t) ≥ θF )l
′ ′V(x , y , t) = 
,



0
(otherwise)(5.12)

where k2 is an integer to set the cut-oﬀ frequency fk2 = f0 k2 /K2 . The peak frequency in
the ROI image lf pk (l = 0,· · ·, L−1) for the l-th labeled honeybee region is determined bysearching the maximum value in the frequency range above the cut-oﬀ frequency fk2 asfollows:
k pk (x′ , y′ , t) = arg max lF K (x′ , y′ , t),lk2 ≤K≤K2′ −1(5.13)

where l k pk (x′ , y′ , t) (l = 0,· · ·, L−1) is an integer to set the peak frequency lf pk = f0 · lk pk /K2
at the pixel (x′ , y′ ) in the ROI image.
To quantify whether there is a certain vibration component of honeybee ﬂapping,
the averaged frequency response lF(x′ , y′ , t) (l = 0,· · ·, L−1) is computed by averaging the
TFRs at the extracted vibration pixels in the ROI images as follows:


1 Xl


V(x′ ,y′, t)lF(x′ ,y′ ,t) (l S (t) ≥ θS )l S (t)l(x′,y′ )∈lO(t)F̄(t) =
,






∅
(otherwise)(5.14)

where l S (t) (l = 0,· · ·, L−1) is the number of the vibration pixels in the l-th labeled honeybee
region l O(t),lS (t) =Xl

V(x′ ,y′, t),(5.15)(x′,y′ )∈lO(t)and θS is a threshold to determine whether the wing-ﬂapping of honeybees is observable.

71

5.4 SENSITIVITY EVALUATION

Independently of the brightness amplitude and the low frequency component in the
background scene, the normalized TFR lA(t) (l = 0,· · ·, L−1) is computed in the range of
their ﬂapping frequencies by normalizing the cut-oﬀ TFR image,l ′F̄k2 (t) = (0, · · · , 0, lF̄k′ 2 (t), · · · , lF̄ K′ ′ −1 ),2(5.16)

whose cut-oﬀ frequency is fk2 , as follows:′′

A(t) = lF̄k2 (t)/ lF̄k2 (t) .

l(5.17)

Assuming that the ﬂapping frequency of a honeybee is larger than a certain frequency fH , the ﬂying activity of the l-th labeled honeybee is computed as its ﬂapping frequency l fA (t) (l = 0,· · ·, L−1) by detecting the peak frequency of lA(t) = (lA0 (t), · · · , lAK2′ −1 (t))
in the following frequency range,lfA (t) =f0
arg max lAK (t),K2 kH ≤K≤K2′ −1(5.18)

where kH is an integer to set the frequency fH = f0 kH /K2 , which determines the frequency
range to be searched.

5.4

Sensitivity evaluation
First, we checked how the spatial resolution of an HFR video aﬀects the sensitiv-

ity in our honeybee activity sensing algorithm when the image intensity at every pixel isconverted to a TFR over a period of time. We captured 8-bit color images of 1024×1024
pixels at 500 fps (τ = 2 ms) with a 0.5-ms-camera exposure when 12–14-mm size honeybees ﬂew around their beehive in a cluttered outdoor location. An HFR camera wasinstalled at a distance of 1.5 m from the beehive. The beehive was observed in the depth

72

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

of ﬁeld (DOF) of the HFR camera using 1024×1024 pixels, and the DOF far limit was set
to inﬁnity with a small lens aperture. A 1024×1024 image corresponded to a 0.6×0.6-marea around the beehive; one pixel corresponded to 0.6 mm. In the evaluation, we used
HFR videos for 0.256 s of 64×64, and 256×256 pixels (M = N = 64 and 256), as wellas the original 1024×1024 HFR video (M = N = 1024), by interleaving the original HFR
video at every 16 and 4 pixel, respectively; the pixel interval in each resolution corresponded to 9.6 and 2.4 mm, respectively. The parameters for step (1) are set as K1 =32, fk0 = 31.25 Hz, fk1 = 62.5 Hz, θF = 400, and θH = 0.5, those for step (2) are set
as m = n = 4, 16, and 64, and θP = 1, 5, and 20 when M = N = 64, 256, and 1024,
respectively, and those for step (3) are set as R = 4, 16, and 64, when M = N = 64, 256,and 1024, respectively, K2 = 128, fk2 = 62.5 Hz, θ′F = 400, θS = 1, and fH = 150 Hz.
These parameters are summarized in Table 5.1.
When the frequency range to be inspected is priorly known for wing-ﬂapping ofhoneybees, K1 , fk0 , fk1 , K2 , fk2 , and fH , which are used to determine the frequency range
and the number of samples in computing TFRs with pixel-level STFTs, can be set as universal values. They are independent of image resolution, image contrast and brightness,
and apparent sizes of honeybees in images; therefore, there is no need to determine themeach time the camera is deployed. The parameters m, n, θP , R, and θS depend on theapparent sizes of honeybees in images, corresponding to camera parameters such as image resolution, magniﬁcation ratio, and camera-to-object distance. They are independent
from time-varying image brightness and contrast in ambient sunlight; therefore, there isno need to adjust these parameters once the camera is deployed at a ﬁxed location. The
parameters θF , θH , and θ′F that mainly check the intensities of the pixel-level TFRs, depend on image contrast and brightness, and background movement. Under time-varying
sunlight condition and background scene, they should be adaptively determined such thathoneybee regions are properly extracted as vibration pixels. However, in this study, we
manually determined these parameters for short-span HFR videos when there was no5.4 SENSITIVITY EVALUATION

73

Table 5.1: Parameters used for the experiments in Subsection 5.4.
64×64
256×2561024×1024
step parameter(M=N=64) (M=N=256) (M=N=1024)
K1
32
fk0
31.25 Hz(1)
fk1
62.5 Hz
θF
400
θH
0.5
m(= n)
4
16
64
(2)
θP
1
5
20
R
4
16
64
K2
128
fk2
62.5 Hz
(3)
′
θF
400
θS
1
fH
150 Hzlarge time-variation in image brightness.
Figure 5.3 shows (a) the input images I(x, y, t), (b) the candidate vibration pixels
V(x, y, t), and (c) the locations of selected ROI images l I(x′ , y′ , t) for honeybees. Figure 5.4 shows the magniﬁed ROI images l I(x′ , y′ , t) for “bee 1” and ”bee 2”, which arelabeled in Figure 5.3. Figure 5.5 shows the averaged TFRs lA(t) in the range from 62.5
to 250 Hz with a frequency resolution of 3.91 Hz, which were computed using the 2nd
pixel-level STFTs for these ROI images to determine the ﬂying activities of the labeledhoneybees. “bee 1” was ﬂying at a certain speed in the direction from right to left, and“bee 2” was hovering in the input images. Figure 5.6 plots (a) the temporal changes in theimage intensities at A(960, 294), B(319, 609), and C(270, 50) in the original input imagesof 1024×1024 pixels and (b) those at A′ (32, 32) and B′ (32, 32) in the tracked ROI imagesof 64×64 pixels for 0.2 s. Here, point A is around “bee 1”, point B is around “bee 2”,
and point C is on a wooden box in the background as illustrated in Figure 5.3(a). Points
A′ and B′ correspond to the centers of the tracked ROI images for “bee 1” and “bee 2”,
respectively, as illustrated in the tracked 64×64 ROI images in Figure 5.4. As shown

74

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

in Figure 5.6, no rapid changes were observed in the image intensity at point C in the
background scene. Only two or three pulses in a short-term duration were observed in theimage intensities at points A and B around the honeybees; they were ﬂying at a certain
speed with their wings ﬂapping at hundreds of Herz. In Figure 5.6(b), periodic waves
with a duration of 0.2 s on the envelope waves, which are drifting with the backgroundand illumination conditions around the honeybees, are observed at points A′ and B′ in the
tracked ROI images; no data was plotted when the tracking was failed. This is becausehoneybee tracking enables the long-time observation of periodic changes in the image
intensities at the pixels around ﬂying honeybees.
In all the cases of 64×64, 256×256, and 1024×1024 input images, candidate vibration pixels were extracted around the ﬂying honeybees with pixel-level digital signal
processing in our algorithm. The honeybees were tracked correctly in the selected ROI
images even when they can be hardly observed in a cluttered scene because their apparent sizes were 1 or 2 pixels in the 64×64 input images. The vibration region of “bee 1”
had a comet-like shape, whose size is larger than the actual size of the honeybee because
the honeybee was ﬂying at a high speed during 64 ms in the STFT computation with
K1 = 32. The honeybee’s velocity was used to determine the direction and length of thecomet shape. The vibration region of “bee 2”, hovering at the same place, had a small
and non-directional shape, compared with that of “bee 1”.
In the averaged TFRs for 256×256 and 1024×1024 input images in Figure 5.5, signiﬁcant frequency components were observed around 175 and 235 Hz, respectively. They
corresponded to the diﬀerent wing-ﬂapping frequencies of “bee 1” and “bee 2”, respectively. Here, the peak frequency detection was conducted as a maximum search process,and it was not very sensitive even in the case of a low signal-to-noise ratio. This was
because the small periodic waves corresponding to the wing ﬂapping of the honeybees
had a distinct frequency component at their wing-ﬂapping frequencies, whereas the TFRs
of the drifting envelope waves were almost ﬂat in the frequency range for peak frequency

5.5 HONEYBEE TRACKING EXPERIMENT

75

detection because they did not have any speciﬁc frequency component. Their peak frequencies were not apparently observed in both the averaged TFRs for “bee 1” and “bee 2”
for 64×64 input images because the locations of the selected ROIs ﬂuctuated pixelwise at
every frame. These ﬂuctuations are not negligible in pixel-level digital signal processing
when one pixel is similar or larger than the size of honeybee in the 64×64 input images.
The spatial resolution of a honeybee in an image was inversely proportional to its distancefrom the camera. Honeybees at a distance of 24 m from the camera could be detected in a
1024×1024 image using pixel-level STFTs, which corresponded to a 9.6×9.6 m area with
a 9.7 mm pixel pitch. This was because honeybees at a distance of 1.5 m from the camera
were detected in 64×64 images, which corresponded to a 0.6×0.6 m area with the same
pixel pitch.
These evaluation results indicate that honeybee tracking can be robustly executed byextracting around ﬂying honeybees in an HFR video using the 1st pixel-level STFTs forvibration source localization, independent of the image resolution. In the meantime, input
images with a certain spatial resolution, in which honeybees are apparently observed, arerequired for the 2nd pixel-level STFTs to stably obtain the frequency responses of ﬂying
honeybees to determine their ﬂying activities.

5.5

Honeybee tracking experiment
To verify the trajectories and ﬂying activities of honeybees in outdoor scene, 1024×1024(M = N = 1024) input images, captured at 500 fps with a 0.5-ms camera exposure for
14 seconds, were analyzed with our activity sensing algorithm. We captured the 500-fps
video of 1024×1024 pixels for the same scene in a 0.6×0.6-m-area around their beehive
in a cluttered outdoor location as that had been analyzed in Subsection 5.4; one pixel
corresponded to 0.6 mm. The parameters in our algorithm were set to the same ones as
those set when M = N = 1024 in Subsection 5.4. In the experiment, most of the ﬂying

76

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

honeybees were wing-ﬂapping at frequencies below the Nyquist frequency at 250 Hz.
Figure 5.7 shows (a) the 1024×1024 input image I(x, y, t) at t = 0 s, (b) the xytrajectories of the labeled honeybees lx(t) for t = 0–14 s, and (c) the ﬂapping frequenciesof the labeled honeybees l fA (t). In (b), the labels of “bee 1”–“bee 10” indicate the indexes in the labeling with moment feature accumulation. These indexes were updated as
those of their nearest honeybees at the previous frame; abandoned indexes that were notassigned at the current frame, were reused when a newly labeled honeybee did not correspond with any labeled honeybees at previous frame. Each xy-trajectory involves many
temporally-disconnected trajectories of diﬀerent honeybees, because the duration times of
their trajectories were mostly less than one second and several trajectories were crossed
over each other. The peak frequencies of the labeled honeybees in (c) were color-mapped
as their ﬂapping frequencies along their trajectories.
It can be observed that there are many complex-shape trajectories of honeybees
around their beehive, whereas there are gentle curved trajectories of honeybees in thehorizontal direction at distant places from their beehive. Here, these trajectories did not
involve the movements of the ﬂying honeybees in the depth direction, because our method
was limited to 2-D measurements in the horizontal and vertical directions. The complexshape trajectories corresponded to those of honeybees hovering around their beehive, andthe ﬂapping frequencies on the trajectories were mostly over 220 Hz; most of them entered
their beehive with hovering for a certain time when they ﬂew back to their beehive. In
this experiment, most of the gentle curve trajectories at distant places from their beehive
corresponded to the honeybees which ﬂew back to their beehive. The ﬂapping frequencies on the horizontal trajectories at distant places were mostly less than 200 Hz, because
the ﬂapping frequencies of honeybees decreased for them to land at their beehives. The
ﬂapping frequencies were uncertain at several segments on the honeybee trajectories dueto the unstable TFRs in brightness with the crosstalk background patterns and the transparency of honeybee wings. We note that the wing-ﬂapping of a honeybee at a frequency

5.5 HONEYBEE TRACKING EXPERIMENT

77

above the Nyquist frequency is undersampled as the low-frequency alias. If a honeybee
is ﬂying with wing-ﬂapping in the range from 250 to 300 Hz in a 500-fps video, it is detected as an undersampled component in the range from 200 to 250 Hz; it is distinctively
separable from the background DC component in vibration source localization.
Figure 5.8 shows (a) the 64×64-input images l (x′ , y′ , t) in the tracked ROI around“bee 2”, (b) the vibration pixels l V(x′ , y′ , t) and (c) the peak frequencies l f pk (x′ , y′ , t) inthe tracked ROI, (d) the time-transient averaged TFR lA(t) in the range of 62.5 to 250 Hz,(e) the apparent ﬂying speed |lu(t)|=|d lx(t)/dt| and ﬂapping frequency lfA (t) of the labeledhoneybee, and (f) the ﬂapping frequency color-mapped on its xy-trajectory when “bee 2”
ﬂew back to its beehive from right to left for t = 6.36–6.96 s; (a)–(c) were obtained atintervals of 0.20 s, and the moving averages of the speeds and ﬂapping frequencies amongcontinuous 10 frames were plotted in (e).
In Figure 5.8, the apparent ﬂying speed of the honeybee decreased from 3.3×103 pixel/s
at t = 6.36 s to 2.7×102 pixel/s at t = 6.96 s in the input images. The apparent speeds
correspond to 2.0 and 0.14 m/s, respectively, around its beehive when the velocity components in the depth direction are ignored; the actual speeds of the honeybees ﬂying in
3-D space are larger than these apparent speeds. It can be observed that a ﬂying honeybee was always in the tracked 64×64 ROI images, and the vibration pixels were detectedaround its ﬂapping wings by detecting the peak frequencies in the tracked ROI. Comparedwith the results in Figure 5.3(b), the vibration pixels have no comet-like shape because
the latency eﬀect in digital signal processing when observing fast-ﬂying honeybees was
remarkably reduced in the tracked ROI images by canceling the translational movementsof honeybees in the original input images. Corresponding to the wing-ﬂapping frequencyof the honeybee, a signiﬁcant frequency component around 175 Hz was observed in the
averaged TFR when the honeybee was ﬂying at 1 m/s or more. Its peak frequency was
ﬂuctuated for t = 6.65–6.85 s because the crosstalk in the averaged TFR was generated
by apparently-moving background patterns in the tracked ROI images.

78

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

In a similar manner, Figure 5.9 and 5.10 show (a) the tracked ROI images, (b) thevibration pixels, (c) the peak frequencies, (d) the averaged TFR, (e) the honeybee speed
and ﬂapping frequency, and (f) the xy-trajectory with ﬂapping frequency when “bee 1”
was hovering around its beehive for t = 7.64–8.24 s, and when “bee 1” took oﬀ and
ﬂew away from its beehive for t = 8.82–9.12 s, respectively; the honeybee labeled as“bee 1” in Figure 5.9 was diﬀerent from that labeled as “bee 1” in Figure 5.10. Whenthe velocity components in the depth direction are ignored, the apparent ﬂying speed ofa honeybee when hovering in Figure 5.9 was less than 1.6×102 pixel/s, corresponding
to 0.10 m/s in the actual scale. The apparent speed when ﬂying away in Figure 5.10
increased from no translation at t = 8.82 s to 7.2×103 pixel/s at t = 9.12 s in input images,
corresponding to 4.3 m/s in the actual scale. In both cases, ﬂying honeybees were stably
tracked in the ROI images, and the vibration pixels were detected around their ﬂappingwings as well as those in Figure 5.8. The ﬂapping frequency when a honeybee washovering was around 235 Hz in Figure 5.9, and that when a honeybee ﬂew away from its
beehive in Figure 5.10 increased from 150 to 205 Hz. The peak frequency in Figure 5.10
was misdetected for t = 8.82–8.85 s and 8.92–8.98 s with background crosstalk in the
tracked ROI images. It can be observed that the ﬂapping frequency of a honeybee timevaried depending on its ﬂight status; the largest ﬂapping frequency was observed when
the honeybee was hovering, followed by that when the honeybee was ﬂying away fromits beehive, and that when the honeybee was ﬂying back to its beehive.
These experimental results show that our HFR-video-based algorithm has pixellevel acute sensitivity to localize and track small ﬂying honeybees in a cluttered outdoor
scene where it is challenging for human eyes to discern them. It has the potential to quantify individually the wing-ﬂapping frequencies of ﬂying honeybees at audio frequency
level, because it inspects their time-varying TFRs in all the tracked ROIs as their ﬂight
activity indexes.

5.6 CONCLUDING REMARKS

5.6

79

Concluding remarks
In this chapter, we proposed a novel honeybee activity sensing algorithm that can

track and quantify the vibration properties in brightness at the pixel around wing-ﬂappinghoneybees by executing two-step pixel-level STFTs for an HFR video. We conducteda tracking experiment for honeybees ﬂying around their beehive in a cluttered scene byanalyzing 1024×1024 images captured at 500 fps, and demonstrated the discriminative
ability of the HFR video to detect the ﬂapping frequencies of honeybees at hundreds ofHertz in the frequency range below 250 Hz as their ﬂight activities, as well as its pixellevel localization ability for individual honeybees. Considering the frequency spectrum
peaks are broad and their time-varying peak frequencies when diﬀerent types of ﬂying
insects are observed in an HFR video, we intend to improve our algorithm to distinguishhoneybees by collecting and analyzing both their appearances and time-varying frequency
spectrums using machine learning technology for pattern recognition. We also intend touse an HFR video at thousands of fps such that the Nyquist frequency is much larger than
the wing-ﬂapping frequencies of insects to be observed, and develop a real-time HFRcamera dedicated to long-time honeybee activity sensing by using graphic processing
units to accelerate the computational speed of pixel-level STFTs in our algorithm.

80

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

宋安宕季宬宱害宸宷季宬宰室宪宨( , , )

= 2

= 3

= 1

weight宦室宱宧宬宧室宷宨
容宬宥宵室宷宬宲宱季害宬宻宨宯家brightness

frequency

×

power

vibration

frequency

害宬宻宨宯孰宯宨容宨宯宖宗安宗家

no vibration

害宬宻宨宯孰宯宨容宨宯
宋宓安家

power

brightness

time

frequencytime宷宬宰宨孰容室宵宼宬宱宪季家宬宪宱室宯季宩宵宨宴宸宨宱宦宼季宵宨家害宲宱家宨季(a) pixel-level vibration source localization
宕宒完季宬宰室宪宨家季宷宲季宥宨季宨宻室宰宬宱宨宧宦宨宯宯孰宥室家宨宧季宯室宥宨宯宬宱宪孴

宷宵室宦宮宬宱宪季宺宬宷宫季
宰宲宰宨宱宷季宩宨室宷宸宵宨
室宦宦宸宰宸宯室宷宬宲宱

孵

孵

孴

孴

孵

孶

孶

孶

⋯

(b) individual honeybee labeling宩宵宨宴宸宨宱宦宼季宵宨家害宲宱家宨季

孴

孴

孴

power

宷宵室宦宮宨宧季宕宒完季宬宰室 宨家 宩宲宵季宫宲宱宨宼宥宨宨季孴( , ,害宬宻宨宯孰宯宨容宨宯季宖宗安宗家季害宬宻宨宯孰宯宨容宨宯季宖宗安宗家季宩宲宵季宕宒完季宬宰室宪宨季孴害宬宻宨宯孰宯宨容宨宯季宖宗安宗家季宩宲宵季宕宒完季宬宰室宪宨季孵
宩宲宵季宕宒完季宬宰室宪宨季孶

⋯

−1 )

宷宵室宭宨宦宷宲宵宬宨家季宲宩季宫宲宱宨宼宥宨宨家孴孶宩宵宨宴宸宨宱宦宼季实宋宽宠孵)宩宵宨宴宸宨宱宦宼季实宋宽宠( , ,宩宵宨宴宸宨宱宦宼季实宋宽宠( , , 1)

frequency宫宲宱宨宼宥宨宨季孴
宼宥

孶孳孳
孵孳孳
孴孳孳
孳

孳

孴

孵

孶

孷

學
宷宬宰宨季实家宠

孹

孳

孴

孵

孶

孷

學
宷宬宰宨季实家宠

孹

孳

孴

孵

孶

孷

學
宷宬宰宨季实家宠

孹

孶孳孳

孺

孻

孼

孴孳

孺

孻

孼

孴孳

孺

孻

孼

孴孳

宫宲宱宨宼宥宨宨季孵

孵孳孳
孴孳孳
孳
孶孳孳

宫宲宱宨宼宥宨宨季孶

孵孳孳
孴孳孳
孳

宫宲宱宨宼宥宨宨季s宼宬宱宪季室宦宷宬容宬宷宬宨家(c) frequency response computation for tracked honeybees
Figure 5.1: Concept of pixel-level ﬂight activity sensing method.

81

5.6 CONCLUDING REMARKS

power

time

frequency

frequency

宩宵宨宴宸宨宱宦宼季宵宨家害宲宱家宨季brightnesstime

家宨容宨宵室宯季宺室容宨家宲宱宨孰宲宵孰宱宲季宺室容宨

time

power

brightness

容宬宥宵室宷宬宲宱

powerbrightness

宬宧宨室宯季宦室家宨

宩宵宨宴宸宨宱宦宼季宵宨家害宲宱家宨季frequency

宩宵宨宴宸宨宱宦宼季宵宨家害宲宱家宨季(a) one or no wave宬宱害宸宷季宬宰室宪宨家害宬宻宨宯孰宯宨容宨宯
宖宗安宗家

宯室宷宨宱宦宼季宨⒐宨宦宷

⋯( , , 1)

( , ,

−1 )

( , ,

)

宪宫宲家宷季容宬宥宵室宷宬宲宱季害宬宻宨宯家(b) latency eﬀect
Figure 5.2: Problems in pixel-level digital signal processing when observing fastﬂying honeybees.

82

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING宆 孫孵孺孳孯學孳孬宄 孫孼孹孳孯孵孼孷孬

宅 孫孶孴孼孯孹孳孼孬孴孳孵孷宻孴孳孵孷季孫宐宀宑宀孴孳孵孷孬孵學孹宻孵學孹季孫宐宀宑宀孵學孹孬孹孷宻孹孷季孫宐宀宑宀孹孷孬(a) input images

孴孳孵孷宻孴孳孵孷季孫宐宀宑宀孴孳孵孷孬孵學孹宻孵學孹季孫宐宀宑宀孵學孹孬孹孷宻孹孷季孫宐宀宑宀孹孷孬(b) vibration pixels

孴孳孵孷宻孴孳孵孷季孫宐宀宑宀孴孳孵孷孬孵學孹宻孵學孹季孫宐宀宑宀孵學孹孬孹孷宻孹孷季孫宐宀宑宀孹孷孬(c) locations of selected ROIs
Figure 5.3: Input images, candidate vibration pixels, and selected ROIs for honeybees in sensitivity evaluation.

83

5.6 CONCLUDING REMARKS宄尋 孫孶孵孯孶孵孬孹孷宻孹孷 宕宒完 孫宐宀宑宀孴孳孵孷孬孴孹宻孴孹 宕宒完 孫宐宀宑宀孵學孹孬孷宻孷 宕宒完 孫宐宀宑宀孹孷孬(a) bee 1

宅尋 孫孶孵孯孶孵孬孹孷宻孹孷 宕宒完 孫宐宀宑宀孴孳孵孷孬孴孹宻孴孹 宕宒完 孫宐宀宑宀孵學孹孬孷宻孷 宕宒完 孫宐宀宑宀孹孷孬(b) bee 2
Figure 5.4: Magniﬁed ROI images in sensitivity evaluation.

84

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

孳孱孳孷

孹孷宻孹孷

宓宲宺宨宵

孳孱孳孶

孵學孹宻孵學孹孴孳孵孷宻孴孳孵孷

孳孱孳孵
孳孱孳孴
孳

孹孵孱學

孻孺孱學

孴孴孵孱學

孴孶孺孱學
孴孹孵孱學
安宵宨宴宸宨宱宦宼季实宋宽宠

孴孻孺孱學

孵孴孵孱學

孵孶孺孱學

孴孻孺孱學

孵孴孵孱學

孵孶孺孱學

(a) bee 1
孳孱孳學

孹孷宻孹孷

宓宲宺宨宵

孳孱孳孷

孵學孹宻孵學孹孴孳孵孷宻孴孳孵孷

孳孱孳孶
孳孱孳孵
孳孱孳孴
孳

孹孵孱學

孻孺孱學

孴孴孵孱學

孴孶孺孱學
孴孹孵孱學
安宵宨宴宸宨宱宦宼季实宋宽宠(b) bee 2
Figure 5.5: Averaged TFRs in the selected ROI images in sensitivity evaluation.

85

5.6 CONCLUDING REMARKS宄季孫孼孹孳孯孵孼孷孬孵學孳

宅季孫孶孴孼孯孹孳孼孬宆季孫孵孺孳孯學孳孬完宱宷宨宱家宬宷宼

孵孳孳
孴學孳
孴孳孳
學孳
孳

孳

孳孱孳孵

孳孱孳孷

孳孱孳孹

孳孱孳孻

孳孱孴

孳孱孴孵

孳孱孴孷

孳孱孴孹

孳孱孴孻

孳孱孵

宗宬宰宨季实家宠

(a) original input images
宄孪季孫孶孵孯孶孵孬孵學孳宅孪季孫孶孵孯孶孵孬完宱宷宨宱家宬宷宼

孵孳孳
孴學孳
孴孳孳
學孳
孳

孳

孳孱孳孵

孳孱孳孷

孳孱孳孹

孳孱孳孻

孳孱孴

孳孱孴孵

孳孱孴孷

孳孱孴孹

孳孱孴孻

孳孱孵

宗宬宰宨季实家宠

(b) tracked ROI images
Figure 5.6: Temporal changes in image intensities in sensitivity evaluation.

Figure 5.7: xy-trajectories and ﬂapping frequencies of labeled honeybees in 14 seconds.

86

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

宷季宀季孹孱孶孹季家

宷季宀季孹孱學孹季家

宷季宀季孹孱孺孹季家

宷季宀季孹孱孼孹季家(a) tracked ROI images

宷季宀季孹孱孶孹季家

宷季宀季孹孱學孹季家

宷季宀季孹孱孺孹季家

宷季宀季孹孱孼孹季家(b) vibration pixels in the tracked ROI
孹孵孱學宷季宀季孹孱孶孹季家

宷季宀季孹孱學孹季家

孵學孳季实宋宽宠宷季宀季孹孱孺孹季家

宷季宀季孹孱孼孹季家(c) peak frequencies in the tracked ROI孳宩宵宨宴宸宨宱宦宼季实宋宽宠

孵學孳

害宲宺宨宵

孴孻孺孱學
孴孵學
孹孵孱學

孹孱孶孹

孹孱孷孹

孹孱學孹

孹孱孹孹
宷宬宰宨季实家宠孹孱孺孹(d) time-transient averaged TFR

孹孱孻孹

宰室宻

室害害室宵宨宱宷季s宼宬宱宪季家害宨宨宧容宨宯宲宦宬宷宼客室宰害宯宬宷宸宧宨

家害宨宨宧季实害宬宻宨宯孲家宠

學孳孳孳
學孳孳孳

害宨室宮季宩宵宨宴宸宨宱宦宼

孵學孳

孷孳孳孳
孷孳孳孳

孵孳孳

孶孳孳孳
孶孳孳孳

孴學孳

孵孳孳孳
孵孳孳孳

孴孳孳

孴孳孳孳
孴孳孳孳
孳孳
孹孱孶孹

孹孱孷孹

孹孱學孹

孹孱孹孹

宷宬宰宨季实家宠

孹孱孺孹

孹孱孻孹

孹孱孼孹

學孳

害宨室宮季宩宵宨宴宸宨宱宦宼季实宋宽宠

87

5.6 CONCLUDING REMARKS(e) apparent ﬂying speed and ﬂapping frequency
孵學孳孵孳孳宷季宀季孹孱孶孹季家
宷季宀季孹孱學孹季家

孴學孳
孴孳孳
學孳

宷季宀季孹孱孺孹季家孳实宋宽宠

宷季宀季孹孱孼孹季家(f) ﬂapping frequency mapped on 0.60-s xy-trajectory
Figure 5.8: Tracked ROI images and ﬂying activities when “bee 2” ﬂew back to its
beehive for t = 6.36–6.96 s.

88

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

宷季宀季孺孱孹孷季家

宷季宀季孺孱孻孷季家

宷季宀季孻孱孳孷 家

宷季宀季孻孱孵孷季家(a) tracked ROI images

宷季宀季孺孱孹孷季家

宷季宀季孺孱孻孷季家

宷季宀季孻孱孳孷 家

宷季宀季孻孱孵孷季家(b) vibration pixels in the tracked ROI
孹孵孱學宷季宀季孺孱孹孷季家

宷季宀季孺孱孻孷季家

孵學孳季实宋宽宠宷季宀季孻孱孳孷 家

宷季宀季孻孱孵孷季家(c) peak frequencies in the tracked ROI孳宩宵宨宴宸宨宱宦宼季实宋宽宠

孵學孳

害宲宺宨宵

孴孻孺孱學
孴孵學
孹孵孱學

孺孱孹孷

孺孱孺孷

孺孱孻孷

孺孱孼孷
宷宬宰宨季实家宠孻孱孳孷(d) time-transient averaged TFR

孻孱孴孷

宰室宻

室害害室宵宨宱宷季s宼宬宱宪季家害宨宨宧容宨宯宲宦宬宷宼客室宰害宯宬宷宸宧宨

家害宨宨宧季实害宬宻宨宯孲家宠

學孳孳孳
學孳孳孳

害宨室宮季宩宵宨宴宸宨宱宦宼

孵學孳

孷孳孳孳
孷孳孳孳

孵孳孳

孶孳孳孳
孶孳孳孳

孴學孳

孵孳孳孳
孵孳孳孳

孴孳孳

孴孳孳孳
孴孳孳孳
孳孳
孺孱孹孷

孺孱孺孷

孺孱孻孷

孺孱孼孷

孻孱孳孷

孻孱孴孷

孻孱孵孷

學孳

害宨室宮季宩宵宨宴宸宨宱宦宼季实宋宽宠

89

5.6 CONCLUDING REMARKS宷宬宰宨季实家宠

(e) apparent ﬂying speed and ﬂapping frequency
孵學孳
孵孳孳
孴學孳
孴孳孳
學孳

宷季宀季孻孱孳孷季家
宷季宀季孻孱孵孷季家

宷季宀季孺孱孹孷季家孳实宋宽宠宷宷季宀季孺孱孻孷季家(f) ﬂapping frequency mapped on 0.60-s xy-trajectory
Figure 5.9: Tracked ROI images and ﬂying activities when “bee 1” was hovering
around its beehive for t = 7.64–8.24 s.

90

CHAPTER 5. HFR-VIDEO-BASED HONEYBEE ACTIVITY SENSING

宷季宀季孻孱孻孵季家

宷季宀季孻孱孼孵季家

宷季宀季孼孱孳孵季家

宷季宀季孼孱孴孵季家(a) tracked ROI images

宷季宀季孻孱孻孵季家

宷季宀季孻孱孼孵季家

宷季宀季孼孱孳孵季家

宷季宀季孼孱孴孵季家(b) vibration pixels in the tracked ROI
孹孵孱學宷季宀季孻孱孻孵季家

宷季宀季孻孱孼孵季家

孵學孳季实宋宽宠宷季宀季孼孱孳孵季家

宷季宀季孼孱孴孵季家(c) peak frequencies in the tracked ROI孳宩宵宨宴宸宨宱宦宼季实宋宽宠

孵學孳

害宲宺宨宵

孴孻孺孱學
孴孵學
孹孵孱學

孻孱孻孵

孻孱孻孺

孻孱孼孵

孻孱孼孺

宷宬宰宨季实家宠

孼孱孳孵(d) time-transient averaged TFR

孼孱孳孺

宰室宻

室害害室宵宨宱宷季s宼宬宱宪季家害宨宨宧容宨宯宲宦宬宷宼客室宰害宯宬宷宸宧宨

家害宨宨宧季实害宬宻宨宯孲家宠

孴孳孳孳孳
孴孳孳孳孳

害宨室宮季宩宵宨宴宸宨宱宦宼

孻孳孳孳
孻孳孳孳

孵孳孳

孹孳孳孳
孹孳孳孳

孴學孳

孷孳孳孳
孷孳孳孳

孴孳孳

孵孳孳孳
孵孳孳孳
孳孳
孻孱孻孵

孵學孳

孻孱孻孺

孻孱孼孵

孻孱孼孺

孼孱孳孵

孼孱孳孺

孼孱孴孵

學孳

害宨室宮季宩宵宨宴宸宨宱宦宼季实宋宽宠

91

5.6 CONCLUDING REMARKS宷宬宰宨季实家宠

(e) apparent ﬂying speed and ﬂapping frequency

宷季宀季孼孱孴孵季家

宷季宀季孻孱孻孵季家
宷季宀季孻孱孼孵季家家

宷季宀季孼孱孳孵季家

孵學孳
孵孳孳
孴學孳
孴孳孳
學孳
孳
实宋宽宠

(f) ﬂapping frequency mapped on 0.30-s xy-trajectory
Figure 5.10: Tracked ROI images and ﬂying activities when “bee 1” ﬂew away fromits beehive for t = 8.82–9.12 s.

Chapter 6
Conclusions
In this study, to realize the visualization of vibration information from target objects
whose dynamics property can be obtained by converting the time-varying signal measured
by an audio or vibration sensor into a frequency in the frequency domain. Moreover, wedeveloped a high-speed-vision-based pixel-level vibration spectrum imaging system inwhich STFTs and digital ﬁlters were implemented in parallel for the time-varying signalsof all the pixels.
We proposed three concepts for the pixel-level vibration spectrum imaging. The
total concept is the pixel-level vibration spectrum imaging method. When a vision system
operates at a suﬃciently high frame rate, temporal periodic changes can be observedin the brightness signals at the pixels around vibrating objects corresponding to theirvibration’s dynamic properties in the audio frequency range. The brightness signals’
dynamic properties can be quantiﬁed in the form of a pixel-level distribution by converting
the brightness signals for all the pixels of the HFR images into their temporal frequency
responses using STFT. This method expands two types of applications for ﬂying objects.
The ﬁrst application proposes a fast algorithm for vision-based vibration source localization for ﬂying multi-copter that can detect vibration sources at hundreds of Hertz byinspecting time-varying brightness signals at each pixel in HFR images. Our algorithm
can signiﬁcantly reduce the computational complexity of pixel-level digital ﬁlters for vi-

93

94

CHAPTER 6. CONCLUSIONS

bration source localization by virtually adjusting the sampling rate to twice the vibration
frequency of a target object to be tracked using downsampled HFR images with frameinterpolation while locking in vibration frequency for tracking objects.
The second application is pixel-level ﬂight activity sensing method by utilizing
pixel-level vibration source localization and vibration spectrum imaging for ﬂying honeybee. This concept shows the ﬂight activity sensing that can track and quantify thevibration properties in brightness at the pixel around wing-ﬂapping honeybees by executing two-step pixel-level STFTs for an HFR video. In several experimental results, ﬂying
honeybees in 1024×1024 images captured at 500 fps demonstrate that, even when the image region of the honeybee is low resolution, robust tracking can be achieved with activity
sensing in the frequency range up to 250 Hz, which is the Nyquist frequency when weanalyze a 500 fps video. The trajectory of a detected and labeled honeybee is monitoredwith the peak frequency of these tracked STFTs as its ﬂight activity.
This pixel-level STFT function for high-frame-rate images at 1000 fps or greater isaccelerated through its parallel-implementation on a GPU-based high-speed vision system, and vibration spectrum images of 512×512 pixels can be obtained in real-time atmore than 50 fps. Its eﬀectiveness was shown through several veriﬁcation experiments,
including the vibration results for a fan, multi-copter, and guitar string, ﬂying honeybee,
which had rotating or vibrating frequencies in the range of hundreds of Hz. These results
can show the possibility of applying for a real ﬁled monitoring system.
The following issues remain to be solved in the future. The detection frequency
resolution is an issue when they use high-speed vision at hundreds and thousands of fps.
As we mentioned the Nyquist frequency and STFT samples, we determined frequency
resolution between these relationships. After detecting target frequency in the proposed
method, we implemented frame-locking into a target frequency to adjust the high accuracy of frequency resolution using the down-sampling method. Although we succeeded
in developing a real-time vibration imaging system at 2000 fps, there is much higher

95

frequency’s phenomena that relate a natural frequency, a resonant frequency, and a chattering frequency and so on than 1000 Hz in several ﬁelds such as milling process andmotor control process. Then, we need to plan to develop a real-time ultrafast vibration
imaging system at 10000 fps or more. We also intend to develop a real-time HFR camera
dedicated to long-time sensing by using GPU to accelerate the computational speed ofpixel-level STFTs in our algorithm. Furthermore, we also plan to research mechanicalvibration analysis, structural inspection, and biological dynamics inspection as appliedresearch on real-time vibration spectrum imaging in future works.

Appendix

Onsite project in TAOYAKA program
Overview of our onsite team projectIn this section, I summarize the activity of the Onsite-Team-Project in the TAOYAKA program, which was held from April 2018 to March 2019. I belong to the TAOYAKA program, which promoted how to learn reverse innovation during my master andPh.D. course at Hiroshima University. In the TAOYAKA program, we made a group of
Onsite-Team-Project consisting of three members, which take up the challenge resolving
multicultural coexistence issues in a disadvantaged area, from a multilateral perspective of
cultural creation, technical creation, and social implementation. Our Onsite-Team-Projectfocused on mobile application based rural development to promote sightseeing in Iinan
town in Shimane prefecture. The sightseeing of hilly and mountainous areas had an essential role in Japan. However, it faces the common issues such as lack of labor and agingof labor and lack of tourists. As the Japanese government would like to implement revitalization eﬀorts outside of Tokyo, this area is drastically declining in population, which
has aﬀected the economy and the social situation. Aging Population of Iinan Town, thepopulation is signiﬁcantly growing in elderly residents. Our team has identiﬁed that Iinan
Town would greatly beneﬁt from a social and economic revitalization project due to its
population decline, low level of tourism, and current economic situation. We proposed97

98

APPENDIX

the social and economic revitalization of Iinan Town community through the development and implementation of a mobile application using Augmented Reality (AR) withembedded video camera and geolocation. It allows players to visit various rural community tourist spots and interact with rural community members by playing game story(currently prototype stage).
This mobile application is explained as follows: When the player launches this application, a mystery map launches the phone’s GPS, which will then provide directions to
the mystery rural town tourist spots. Players can see these cats only by looking throughthis application interface. Diﬀerent interactions allow you to befriend the virtual charactercats. When the player visits the tourist spot and scans the QR code at the tourist location,a cat then appears through AR. The Players visit the tourist spots to learn more about how
to befriend the cats and ﬁnd them by talking with community members. The objective of
this application is for players from all over Japan to visit the rural community to interact
with locals while visiting tourist spots and spend money within the community. Moreover, we propose the sustainable application development through the mobile application
to focus on tourism regarding young generation education and I-turn and U-turn to increase the population. Increasing tourism in rural communities not only boosts the social
interaction among the elderly in the community but will possibly boost the economy of
the community.
For my part as a member of the Onsite-Team-Project, I assisted with the development of a mobile application game using AR technology with integrated geo-location
tagging software to promote tourism, economic growth, and stimulate social interaction
for Iinan Town. Our team cooperated with the University of Texas at Dallas, Narrative
Systems laboratory in the USA.
Through this activity, I suggested that we promote a programming school as thismobile application can be a sustainable development for the local community.
If local people continue to develop this application, they can get the beneﬁts to pro-

99

APPENDIX

Figure A.1: Our prototype mobile application and future extension
mote tourism and expand IT Start-up companies and IT education development throughprogramming schools to stimulate their local economy. When they recruit a game froma local school story after getting programming skills, they can implement a virtual entertainment park with plenty of natural resources.
We purchased embedded systems for education kits and AI sticks such as NVIDIA
jetson tx2 and Intel Movidius stick in order to promote education for programming skills.
We applied YOLO (You Only Look Once) [150], which is a kind of AI technology-based
object recognition and localization. We can easily assign the AR object with this technology because it includes appearance feature-based recognition and localization. Figure 6
shows that our proposed prototype application and future extension of an example ofobject recognition and gender and age classiﬁcation [151] using convolutional neural networks with embedded systems. Local people can expand AR mobile applications withthis equipment because programming skills implemented these technologies. Programming skills can expand the future application.

100

APPENDIX

Future onsite implementation
Through the above activity, we have plans for future onsite implementation combining this study’s technology. Advanced mobile phone camera technology, which cancapture images in hundreds of frames per second, can be used in disadvantaged areas suchas rural communities with none or insuﬃcient electric infrastructures. For example, rural
areas have old structures such as houses and bridges to inspect manually for foundation
stability. Although there are not many infrastructural engineers within these rural areas.In
this situation, having a more straightforward method to inspect the stability of their foundations, infrastructure, and monitoring systems are very important for rural community
members. Introducing an HFR video analysis that can extract pixel-level temporal frequency responses at hundreds of Hz or more, for structure, health inspection, and evenmonitoring, it would be a huge beneﬁt to rural communities in Japan. It indicates thata cellular phone with an HFR video camera can function as a convenient mobile structure inspection system by implementing a cloud service-based software application for
pixel-level temporary frequency analysis via the internet.
As we proposed a mobile application with geolocation, we can distribute eﬃcient
inspection through sightseeing activity using location map in some event. Additionally
I would like to promote developing an automatic honeybee monitoring system using our
proposed method in a disadvantaged area that has plenty of ﬂower resources.
As we explained about CCD in chapter 5.1, CCD has become a severe social problem that causes signiﬁcant economic losses to agricultural business because honeybees
pollinate many crops. The number of honeybees is remarkably decreasing due to various issues such as enemy invasion and some disease, but the demand for honeybee isincreasing in various ﬁelds such as the Agricultural ﬁeld and medical ﬁeld.
Although there is a signiﬁcant demand for honeybee for management, it’s challenging to monitor their activity and condition remotely. The younger generation and elderly
people can also manage honeybee monitoring using automated honeybee activity moni-

101

APPENDIX

Figure A.2: Future onsite implementation for honeybee monitoring
toring technology. Figure 6 shows that our proposed system for honeybee monitoring caninform beekeepers using the mobile application when some incidents happen. When they
can monitor which beehive is active or non-active and which beehive is attacked or not,they can manage the beehive more easily. If they utilize the honeybee monitoring system,
they will make nutritious honey and their agricultural, economic, and health beneﬁts more
eﬃcient. These plans have a sustainable possibility for onsite reverse innovation.

Bibliography[1] Y. Lei, J. Lin, M.J. Zuo, and Z. He, “Condition monitoring and fault diagnosis ofplanetary gearboxes: A review,” Measurement, vol. 48, pp. 292–305, 2014.
[2] C. Li, R. V. Sanchez, G. Zurita, M. Cerrada, and D. Cabrera, “Fault diagnosis forrotating machinery using vibration measurement deep statistical feature learning,”
Sensors, vol. 16, no. 6, pp. 895, 2016.
[3] M. Saﬁzadeh, and S. Latiﬁ, “Using multi-sensor data fusion for vibration fault diagnosis of rolling element bearings by accelerometer and load cell,” Inform. Fusion,
vol. 18, pp. 1–8, 2014.
[4] A. Stanbridge, M. Martarelli, and D. Ewins, “Measuring area vibration mode shapeswith acontinuous-scan LDV,” Measurement, vol. 35, pp. 181–189, 2004.
[5] J. Torano, R. Rodriguez, I. Diego, J. M. Rivas, and M.D. Casal, “FEM models including randomness and its application to the blasting vibrations prediction,” Computers and Geotechnics, vol. 33, no. 1, pp. 15–28, 2006.
[6] X. W. Ye, Y. Q. Ni, T. T. Wai, K. Y. Wong, X. M. Zhang, and F. Xu, “A vision-based
system for dynamic displacement measurement of long-span bridges: algorithm andveriﬁcation,” Smart Struct. Syst., vol. 12, no. 3, pp. 363–379, 2013.
[7] D. Feng, and M. Q. Feng, “Vision-based multipoint displacement measurement forstructural health monitoring,” Struct. Control Health Monit., vol. 23, no. 5, pp. 876–
890, 2015.
103

104

BIBLIOGRAPHY[8] E. Caetano, S. Silva, and J. Bateira, “A vision system for vibration monitoring ofcivil engineering structures,” Exp. Tech., vol. 35, no. 4, pp. 74–82, 2011.
[9] G. Busca, A. Cigada, P. Mazzoleni, E. Zappa, “Vibration monitoring of multiple
bridge points by means of a unique vision-based measuring system,” Exp. Mech.,vol. 54, pp. 255–271, 2014.
[10] X. Dai, X. Shao, Z. Geng, F. Yang, Y. Jiang, and X. He, “Vibration measurement
based on electronic speckle pattern interferometry and radial basis function,” Opt.
Commun., vol. 355, pp. 33–43, 2015.
[11] I. Ishii, I. Ohara, T. Tatebe, and T. Takaki, “1000 fps target tracking using vibrationbased image features,” Proc. IEEE Int. Conf. Robot. Autom., pp. 1837–1842, 2011.
[12] M. Jiang, T. Aoyama, T. Takaki, and I. Ishii. “Pixel-level and robust vibration source
sensing in high-frame-rate video analysis,” Sensors, vol. 16, no.11, pp. 1842, 2016.
[13] M. Jiang, Q. Gu, T. Aoyama, T. Takaki, and I. Ishii. “Real-time vibration source
tracking using high-speed vision,” IEEE Sensors J., vol. 17, pp. 1513–1527, 2017.
[14] I. Ishii, T. Taniguchi, R. Sukenobe, and K. Yamamoto, “Development of high-speed
and real-time vision platform, H3 Vision,” IEEE/RSJ Int. Conf. Intell. Robots andSyst., pp. 3671–3678, 2009.
[15] I. Ishii, T. Tatebe, Q. Gu, Y. Moriue, T. Takaki, and K. Tajima, “2000 fps real-time
vision system with high-frame-rate video recording,” Proc. IEEE Int. Conf. Robot.
Automat., pp. 1536–1541, 2010.
[16] T. Yamazaki, H. Katayama, S. Uehara, A. Nose, M. Kobayashi, S. Shida, M. Odahara, K. Takamiya, Y. Hisamatsu, S. Matsumoto, L. Miyashita, Y. Watanabe, T.
Izawa, Y. Muramatsu, and M. Ishikawa, “A 1 ms high-speed vision chip with 3D-

BIBLIOGRAPHY

105

stacked 140GOPS column-parallel PEs for spatio-temporal image processing,” Abst.
Int. Solid-State Circuits Conf., pp. 82–83, 2017.
[17] A. Sharma, K. Shimasaki, Q. Gu, J. Chen, T. Aoyama, T. Takaki, and I. Ishii, “Super
high-speed vision platform that can process 1024×1024 images in real time at 12500
fps,” Proc. IEEE/SICE Int. Symp. Syst. Integr., pp. 544–549, 2016.
[18] D. Zhang, J. Guo, X. Lei, and C. Zhu, “A High-speed vision-based sensor for dynamic vibration analysis using fast motion extraction algorithms,” Sensors, pp. 16–
572, 2016.
[19] Q. Gu, T. Takaki, and I. Ishii, “A fast multi-object extraction algorithm based oncell-based connected components labeling,” IEICE Trans. Inform. Syst., vol. E95-D,
no. 2, pp. 636–645, 2012.
[20] Q. Gu, T. Takaki, and I. Ishii, “Fast FPGA-based multi-object feature extraction,”
IEEE Trans. Circ. Sys. Video Tech., vol. 23, no. 1, pp. 30–45, Jan. 2013.
[21] K. Okumura, K. Yokoyama, H. Oku, and M. Ishikawa, “1ms auto pan-tilt – video
shooting technology for objects in motion based on saccade mirror with backgroundsubtraction,” Adv. Robot., vol. 29, no. 7, pp. 457–468, 2015.
[22] A. Namiki, K. Hashimoto, and M. Ishikawa, “A hierarchical control architecture forhigh-speed visual servoing,” Int. J. Robot. Res., vol. 22, no. 10–11, pp. 873–888,
2003.
[23] T. Senoo, A. Namiki, and M. Ishikawa, “Ball control in high-speed batting motion
using hybrid trajectory generator,” Proc. IEEE Int. Conf. Robot. Automat., pp. 1762–
1767, 2006.
[24] A. Namiki and N. Ito, “Ball catching in Kendama game by estimating grasp condi-

106

BIBLIOGRAPHY

tions based on a high-speed vision system and tactile sensors,” Proc. IEEE-RAS Int.
Conf. Humanoid Robots, pp. 634–639, 2014.
[25] T. Aoyama, T. Takaki, T. Miura, Q. Gu, and I. Ishii, “Realization of ﬂower stickrotation using robotic arm,” Proc. IEEE/RSJ Int. Conf. Intell. Robots and Syst., pp.
5648–5653, 2015.
[26] L. Wang, R. X. Gao, “Condition monitoring and control for intelligent manufacturing,” Springer, Berlin, 2006.
[27] E. Jantunen, “A summary of methods applied to tool condition monitoring in
drilling,” Int. J. Mach. Tools Manuf., vol. 42, no. 9, pp. 997–1010, 2002.
[28] A. Z. Khan, A. B. Stanbridge, and D.J. Ewins, “Detecting damage in vibrating structures with a scanning LDV,” Opt. Lasers Eng., vol. 32, no. 6, pp. 583–592, 1999.
[29] D. Goyal, B. S. Pabla, “The Vibration Monitoring Methods and Signal Processing
Techniques for Structural Health Monitoring: A Review,” Arch. Computat. Methods
in Eng. vol. 23, pp. 585–594, 2016.
[30] M. Dowell, and G. Sylvester, “Turbomachinery prognostics and health management
via eddy current sensing: current developments,” Proc. IEEE Aerospace Conf., vol3, pp 1–9, 1999.
[31] L. Maurice and M. L. Adams, “Rotating machinery vibration: from analysis to troubleshooting,” CRC Press, 2000.
[32] C. P. Lawson and P. C. Ivey, “Tubomachinery blade vibration amplitude measurement through tip timing with capacitance tip clearance probes,” Sens. Actuators, vol.
118, no. 1, pp. 14–24, 2005.
[33] R. G. Silva, R. L. Reuben, K. J. Baker, and S. J. Wilcox, “Tool wear monitoring ofturning operations by neural network and expert system classiﬁcation of a feature

BIBLIOGRAPHY

107

set generated from multiple sensors,” Mech. Syst. Signal. Process., vol. 12, no. 2,
pp. 319–332, 1998.
[34] Dimla DE, Lister PM “On-line metal cutting tool condition monitoring. I: force andvibration analyses.” Int. J. Mach. Tools Manuf., 2000
[35] Zhou, H., Yuan, Y. and Shi, C., “Object tracking using SIFT features and mean
shift,” Comput. Vision and Image Understanding, vol. 113, no. 3, pp. 345–352, 2009.
[36] S. Pavani, A. Delgado-Gomez, and A. Frangi, “Gaussian weak classiﬁers based on
co-occurring haar-like features for face detection,” Pattern Anal. and Appl., vol. 17,
no. 2, pp. 431–439, 2014.
[37] B. Wu, C. Kao, C. Jen, Y. Li, Y. Chen, and J. Juang, “A relative-discriminativehistogram-of-oriented-gradients-based particle ﬁlter approach to vehicle occlusion
handling and tracking,” IEEE Trans. Ind. Electronics, vol. 61, no. 8, pp. 4228–4237,
2014.
[38] A. Satpathy, X. Jiang, and H. L. Eng, “LBP-based edge-texture features for object
recognition,” IEEE Trans. Image Process., vol. 23, no. 5, pp. 1953–1964, 2014.
[39] S. Sivaraman, and M. M. Trivedi, “Looking at vehicles on the road: A survey ofvision-based vehicle detection, tracking and behavior analysis,” IEEE Trans. Intell.
Transp. Syst., vol. 14, no. 4, pp. 1773–1795, 2013.
[40] S. Rautaray, and A. Agrawal, “Vision based hand gesture recognition for human
computer interaction: a survey,” Artif. Intell. Rev., vol. 43, no. 1, pp. 1–54, 2015.
[41] I. Ishii, T. Taniguchi, K. Yamamoto, and T. Takaki, “High-frame-rate optical ﬂow
system,” IEEE Trans. Circuits Syst. Video Technol., vol. 22, no. 1, pp. 105–112,
2012.

108

BIBLIOGRAPHY[42] I. Ishii, T. Tatebe, Q. Gu, and T. Takaki, “Color-histogram-based tracking at 2000
fps,” J. Electron. Imaging, vol. 21, no. 1, pp. 013010, 2012.
[43] Q. Gu, S. Raut, K. Okumura, T. Aoyama, T. Takaki, and I. Ishii, “Real-time image
mosaicing system using a high-frame-rate video sequence,” J. Robot. Mechatronics,
vol. 27, no. 1, pp. 12–23, 2015.
[44] I. Ishii, T. Ichida, Q. Gu, and T. Takaki. “500-fps face tracking system,” J. Real-time
Image Process., vol. 8, no. 4, pp. 379–388, 2013.
[45] S. Hertega, and J. Liljencrantz, “Measurement of human vocal fold vibrations withlaser triangulation,” Opt. Eng., vol. 40, no. 9, pp.2041–2044, 2001.
[46] J. H. Wu, R. S. Chang, and J. A. Jiang, “A novel pulse measurement system by using
laser triangulation and a CMOS image sensor,” Sensors, vol. 7, no. 12, pp.3366–
3385, 2007.
[47] P. Castellini, M. Martarelli, and E. P. Tomasini, “Laser doppler vibrometry: Development of advanced solutions answering to technology’s needs,” Mech. System andSignal Process., vol. 20, no. 6, pp. 1265–1285, 2006.
[48] H. Oku, I. Ishii, and M. Ishikawa, “Tracking a protozoon using high-speed visualfeedback,” Proc. Int. IEEE-EMBS Conf. Microtechnologies in Medicine and Biology, pp. 156–159, 2000.
[49] Q. Gu, T. Kawahara, T. Aoyama, T. Takaki, I. Ishii, A. Takemoto, and N. Sakamoto,“LOC-based high-throughput cell morphology analysis system,” IEEE Trans. Automat. Sci. Eng., vol. 12, no. 4, pp. 1346–1356, 2015.
[50] H. Yang, Q. Gu, T. Aoyama, T. Takaki, and I. Ishii, “Dynamics-based stereo visualinspection using multidimensional modal analysis,” IEEE Sensors Journal, vol. 13,
no. 12, pp. 4831–4843, 2013.

BIBLIOGRAPHY

109[51] S. Argentieri, P. Danes, and P. Soueres, “A survey on sound source localization in
robotics: Binaural to array processing methods,” Comput. Speech Lang., vol. 34, no.
1, pp. 87–112, 2015.
[52] Y. C. Lu and M. Cooke, “Binaural estimation of sound source distance via thedirect-to-reverberant energy ratio for static and moving sources,” IEEE Trans. Audio, Speech, Lang. Process., vol. 18, no. 7, pp. 1793–1805, 2010.
[53] U. H. Kim, K. Nakadai, and H. G. Okuno, “Improved sound source localization inhorizontal plane for binaural robot audition,” Appl. Intell., vol. 42, no. 1, pp. 63–74,
2015.
[54] X. Zhong, W. Yost, and L. Sun, “Dynamic binaural sound source localization with
ITD cues: Human listeners,” J. Acoust. Soc. Amer., vol. 137, no. 4, pp. 2376–2376,
2015.
[55] M. Raspaud, H. Viste, and G. Evangelista, “Binaural source localization by joint
estimation of ILD and ITD,” IEEE Trans. Audio, Speech, Lang. Process., vol. 18,
no. 1, pp. 68–77, 2010.
[56] C. H. Knapp and G. C. Carter, “The generalized correlation method for estimationof time delay,” IEEE Trans. Acoust. Speech Signal Process., vol. 24, no. 4, pp. 320–
327, 1976.
[57] J. Benesty, “Adaptive eigenvalue decomposition algorithm for passive acoustic
source localization,” J. Acoust. Soc. Amer., vol. 107, no. 1, pp. 384–391, 2000.
[58] J. Chen, J. Benesty, and Y. Huang, “Time delay estimation in room acoustic environments: an overview,” EURASIP J. Appl. Signal Process., vol. 2006, pp. 1–19,
2006.

110

BIBLIOGRAPHY[59] R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE
Trans. Antennas Propag., vol. 34, no. 3, pp. 276–280, 1986.
[60] K. Okutani, T. Yoshida, K. Nakamura, and K. Nakadai, “Outdoor auditory scene
analysis using a moving microphone array embedded in a quadrocopter,” Proc.
IEEE/RSJ Int. Conf. Intell. Robots Syst., pp. 3288–3293, 2012.
[61] L. J. Griﬃths and C. W. Jim, “An alternative approach to linearly constrained adaptive beamforming,” IEEE Trans. Antennas Propag., vol. 30, no. 1, pp. 27–34, 1982.
[62] J. M. Valin, F. Michaud, and J. Rouat, “Robust localization and tracking of simultaneous moving sound sources using beamforming and particle ﬁltering,” Robot.
Auton. Syst., vol. 55, no. 3, pp. 216–228, 2007.
[63] F. Ribeiro, D. Florencio, D. Ba, and C. Zhang, “Geometrically constrained roommodeling with compact microphone arrays,” IEEE Trans. Audio Speech Lang. Process., vol. 20, no. 5, pp. 1449–1460, 2012.
[64] F. Keyrouz, “Advanced binaural sound localization in 3-D for humanoid robots,”
IEEE Trans. Instrum. Meas., vol. 63, no. 9, pp. 2098–2107, 2014.
[65] R. Li, and D. He, “Rotational machine health monitoring and fault detection using
EMD-based acoustic emission feature quantiﬁcation,” IEEE Trans. Instrum. Meas.,
vol. 61, no. 4, pp. 990–1001, 2012.
[66] J. A. Ballesteros, E. Sarradj, M. D. Fernandez, T. Geyer, and M. J. Ballesteros,“Noise source identiﬁcation with beamforming in the pass-by of a car,” Appl.
Acoust., vol. 93, pp. 106–119, 2015.
[67] P. Castellini, A. Sassaroli, A. Paonessa, A. Peiﬀer, and A. Roeder, “Average beamforming in reverberant ﬁelds: Application on helicopter and airplane cockpits,”
Appl. Acoust., vol. 74, no. 1, pp. 198–210, 2013.

111

BIBLIOGRAPHY[68] K. Nakadai, T. Takahashi, H. G. Okuno, H. Nakajima, Y. Hasegawa, and H. Tsujino, “Design and implementatiaon of robot audition system HARK - open source
software for listening to three simultaneous speakers,” Adv. Robot., vol. 24, no. 5-6,
pp. 739–761, 2010.
[69] D. Yang, Z. Wang, B. Li, and X. Lian, “Development and calibration of acoustic
video camera system for moving vehicles,” J. Sound and Vib., vol. 330, no. 11,
pp. 2457–2469, 2011.
[70] B. Barsikow, W. F. King, and E. Pﬁzenmaier, “Wheel/rail noise generated by ahigh-speed train investigated with a line array of microphones,” J. Sound and Vib.,
vol. 118, no. 1, pp. 99–122, 1987.
[71] Y. Li, M. Smith, and X. Zhang, “Measurement and control of aircraft landing gear
broadband noise,” Aerosp. Sci. Technol., vol. 23, no. 1, pp. 213–223, 2012.
[72] T. Pham, and N. Srour, “TTCP AG-6: acoustic detection and tracking of UAVs,”
Defense and Secur., pp. 24–30, 2004.
[73] J. Busset, F. Perrodin, P. Wellig, B. Ott, K. Heutschi, T. Ruhl, and T. Nussbaumer,“Detection and tracking of drones using advanced acoustic cameras,” SPIE Secur. +Defence, pp. 96470F–96470F, 2015.
[74] Dedrone,
USA,“Multi-sensor(2020).

drone

Accessed:June.

warning15,

system.”
2020.

San

francisco,[Online].

Available:

http://www.dedrone.com/en/dronetracker/drone-detection-hardware.
[75] Droneshield, “How droneshield works.” Sydney, Australia, (2020). Accessed: June.
15, 2020. [Online]. Available:
works1.

https://www.droneshield.com/how-droneshield-

112

BIBLIOGRAPHY[76] K. Shimasaki, T. Okamura, M. Jiang, T. Takaki, I. Ishii, and K. Yamamoto, “HFRvideo-based image pattern recognition using pixel-level temporal frequency response matching,” Proc. IEEE Int. Conf. Automat. Sci. Eng., pp. 451–456, 2018.
[77] J. Seo, H. Ma, and T. K. Saha, “A joint vibration and arcing measurement systemfor online condition monitoring of on-load tap changer of power transformer,” IEEE
Trans. Power Del., vol. 32, no. 2, pp. 1031–1038, 2017.
[78] Y. Park, M. Jeong, S. B. Lee, J. A. Antonino-Daviu, and M. Teska, “Inﬂuence of
blade pass frequency vibrations on MCSA-based rotor fault detection of induction
motors,” IEEE T. Ind. Appl., vol. 53, no. 3, pp. 2049–2058, 2017.
[79] W. Pan, X. Kong, and J. Xu, “Measurement and analysis system of vibration forthe detection of insect acoustic signals,” 2016 Asia-Paciﬁc Int. Symp. Electromagn.
Compat. (APEMC), 2016.
[80] N. S. Marne, M. S. Nagmode, and R. D. Komati, “Vibration measurement systemwith accelerometer sensor based on ARM,” Int. J. Emerg. Technol. Adv. Eng., vol.
4, no. 4, 2014.
[81] A. Devillez and D. Daniel, “Tool vibration detection with eddy current sensors inmachining process and computation of stability lobes using fuzzy classiﬁers, Mech.
Sys. Signal Pr., vol. 21, no. 1, pp. 441–456, 2007.
[82] Y. R. Garcia, J. M. Corres, and J. Goicoechea, “Vibration detection using opticalﬁber sensors,” J. Sens., vol. 2010, pp. 1–12 2010.
[83] F. Garoi, P. C. Logofatu, D. Apostol, C. Udrea, and P. Schiopu, “InterferometricVibration displacement measurement,” Physics, vol. 62, pp. 671–677, 2010.
[84] D. Mas, J. Espinosa, J. Perez, C. Illueca, B. Ferrer, and A. B. Roig, “High speed

BIBLIOGRAPHY

113

imaging and algorithms for non-invasive vibrations measurement,” Proc. EVACES
2011 Exp. Vib. Anal. Civ. Eng. Struct., pp. 531–538, 2011.
[85] E. Carrera and A. Pagani, “Free vibration analysis of civil engineering structures bycomponent-wise models,” J. Sound Vib., vol. 333, no. 19, pp. 4597–4620, 2014.
[86] S. Wan, D. Li, J. Gao, R. Roy, and Y. Tong, “Process and knowledge management ina collaborative maintenance planning system for high value machine tools,” Comput.
Ind., vol. 84, pp. 14–24, 2017.
[87] Y. Nie, I. Ishii, K. Yamamoto, K. Orito, and H. Matsuda, “Real-time scratching
behavior quantiﬁcation system for laboratory mice using high-speed vision,” J. RealTime Image Process., vol. 4, pp. 181–190, 2009.
[88] S. Sakuma, K. Kuroda, C. Tsai, W. Fukui, F. Arai, and M. Kaneko, “Red blood cell
fatigue evaluation based on the close-encountering point between extensibility andrecoverability,” Lab on a Chip, vol. 14, pp. 1135–1141, 2014.
[89] Q. Gu, T. Aoyama, T. Takaki, and I. Ishii, “Simultaneous vision-based shape andmotion analysis of cells fast-ﬂowing in a microchannel,” IEEE Trans. Autom. Sci.
Eng., vol. 12, pp. 204–215, 2015.
[90] H. Yang, T. Takaki, and I. Ishii, “Simultaneous Dynamics-Based Visual Inspection
Using Modal Parameter Estimation,” J. Robot. Mech., vol. 23, no. 1, 2011.
[91] T. Aoyama, L. Li, M. Jiang, K. Inoue, T. Takaki, I. Ishii, H. Yang, C. Umemoto,
H. Matsuda, M. Chikaraishi, and A. Fujiwara, “Vibration sensing of a bridge modelusing a multithread active vision system,” IEEE/ASME Trans. Mechatronics, vol. 23,
no 1, pp. 179–189, 2017.
[92] “Anti-drone market by technology (laser, kinetic, and electronics), application (de-114

BIBLIOGRAPHY

tection and detection & disruption), vertical (military & defense, homeland security,and commercial), and geography - global forecast to 2023,” 2017,[93] A. Yilmaz, O. Javed, and M. Shah.“Object tracking: a survey,” ACM Comput. Surv.,
vol. 38, no. 4, Jan. 2006.
[94] A. Rozantsev, V. Lepetit, and P. Fua, “Detecting ﬂying objects using a single moving
camera,” IEEE Trans. Pattern Anal. and Mach. Intell., vol. 39, no. 5, pp. 879–892,
2017.
[95] R. Chandra, R., Menon, L., Dagum, D., Kohr, D., Maydan, and J. McDonald, “Parallel programming in openMP,” Morgan Kaufmann Publishers, San Francisco, CA.,
2000.
[96] G. B. Patricio-Roberto and M. J. O. Campos, “Aspects of landscape and pollinatorswhat is important to bee conservation?,” Diversity, vol. 6, pp. 158–175, 2014.
[97] G. R. Williams, D. R. Tarry, D. van Engelsdorp, M. P. Chauzat, D. L. Cox-Foster,K. S. Delaplane, P. Neumann, J. S. Pettis, R. E. L. Rogers, and D. Shutler, “Colony
collapse disorder in context,” BioEssays, vol. 32, no. 10, pp. 845–846, 2010.
[98] D. van Engelsdorp, J. D. Evans, C. Saegerman, C. Mullin, E. Haubruge, B. K.
Nguyen, M. Frazier, J. Frazier, D. Cox-Foster, Y. Chen, R. Underwood, D. R. Tarpy,
and J. S. Pettis, “Colony collapse disorder: A descriptive study,” PLoS One, vol. 4,
no. 8, pp. 1–17, 2009.
[99] D. van Engelsdorp, D. Caron, J. Hayes, R. Underwood, M. Henson, A. Spleen, M.
Andree, R. Snyder, K. Lee, K. Roccasecca, M. Wilson, J. Wilkes, E. Lengerich, and
J. Pettis “A national survey of managed honey bee 2010-11 winter colony losses
in the USA: results from the bee informed partnership,” J. Apic. Res., vol. 51, pp.
115–124, 2012.

BIBLIOGRAPHY

115[100] D. Goulson, E. Nicholls, C. Botias, and E. L. Rotheray, “Bee declines driven by
combined stress from parasites, pesticides and lack of ﬂowers,” Science, vol. 347,
no. 6229, pp. 1255957, 2015.
[101] P. Neumann and N. L. Carreck, “Honey bee colony losses,” J. Apic. Res., vol. 49,
pp. 1–6, 2010.
[102] B. N. Gates, ”The temperature of the bee colony,” Bull. US Dept. Agric., vol. 96,
pp. 1–29, 1914.
[103] W. G. Meikle and N. Holst, “Application of continuous monitoring of honeybeecolonies,” Apidologie, vol. 46, no. 1, pp. 10–22, 2015.
[104] W. G. Meikle, N. Holst, G. Mercadier, F. Derouane, and R. R. James, “Usingbalances linked to dataloggers to monitor honeybee colonies,” J. Apic. Res., vol. 45,
no. 1, pp. 39–41, 2006.
[105] W. G. Meikle, B. G. Rector, G. Mercadier, and N. Holst, “Within-day variation incontinuous hive weight data as a measure of honey bee colony activity,” Apidologie,
vol. 39, no. 6, pp. 694–707, Dec. 2008.
[106] S. C. Thoenes and S. L. Buchmann, “Colony abandonment by adult honey bees:
A behavioral response to high tracheal mite infestation?,” J. Apic. Res. vol. 31, no.
3–4, pp. 167–168, 1992.
[107] E. E. Southwick and J. N. Mugaas, “A hypothetical homeotherm: The honey bee
hive,” Comp. Biochem. Physiol. vol. 40, no. 4, pp. 935–944, 1971.
[108] K. V. Nerum and H. Buelens, “Hypoxiacontrolled winter metabolism in honeybees(Apis mellifera),” Comp. Biochem. Physiol. vol. 117, no. 4, pp. 445–455, 1997.
[109] F. Kronenberg and H. C. Heller, “Colonial thermoregulation in honey bees (Apismellifera),” J. Comp. Physiol. vol. 148, no. 1, pp. 65–76, 1982.

116

BIBLIOGRAPHY[110] H. Human, S. W. Nicolson, and V. Dietemann, “Do honeybees,” Apis mellifera
scutellata, vol. 93, no. 8, pp. 397–401, 2006.
[111] K. V. Frisch. “The dance language and orientation of bees,” Harvard University
Press, Cambridge, Massachusetts, 1967.
[112] L. Fahrenholz, I. Lamprecht, and B. Schricker, “Calorimetric investigations of thediﬀerent castes of honey bees, apis mellifera carnica,” J. Comp. Physiol. B, vol. 162,
no. 2, pp. 119–130, 1992.
[113] H. A. Toufailia, M. Couvillon, F. Ratnieks, and C. Grüter, “Honey bee waggle
dance communication: Signal meaning and signal noise aﬀect dance follower behaviour,” Behav. Ecol. Sociobiol., vol. 67, no. 4, pp. 549–556, 2013.
[114] M. Hrncir, F. G. Barth, and J. Tautz, “Vibratory and airborne-sound signals inbee communication (Hymenoptera),” Insect Sound Commun. Physiol. Behav. Ecol.
Evol., pp. 421, 2005.
[115] J. C. Nieh and J. Tautz, “Behaviour-locked signal analysis reveals weak 200-300
Hz comb vibrations during the honey bee waggle dance,” J. Exp. Biol. vol. 203, no.
10, pp. 1573–1579, 2000.
[116] M. Bencsik, J. Bencsik, M. Baxter, A, Lucian, J. Romieu, and M. Millet, “Identiﬁcation of the honey bee swarming process by analysing the time course of hive
vibrations,” Comput. Electron. Agric. vol. 76, no. 1, pp. 44–50, 2011.
[117] S. Ferrari, M. Silva, M., Guarino, and D. Berckmans, “Monitoring of swarming
sounds in bee hives for early detection of the swarming period,” Comput. Electron.
Agric. vol. 64, no. 1, pp. 72–77, 2008.
[118] D. A. Mezquida and J. L. Martinez, “Platform for bee-hives monitoring based on

BIBLIOGRAPHY

117

sound analysis. A perpetual warehouse for swarm’s daily activity span,” J. Agric.
Res., vol. 7, no. 4, pp. 824–828, 2009.
[119] A. R, McLellan, “Honey bee colony weight as an index of honey production andnectar ﬂow: a critical evaluation,” J. Appl. Ecol., vol. 14, no. 2, pp. 401–408, 1997.
[120] M. H. Pham-Delegue, A. Decourtye, L. Kaiser, and J. Devillers, “Behavioural
methods to assess the eﬀects of pesticides on honey bees,” Apidologie vol. 33, no.
5, pp. 425–432, 2002.
[121] A. C. Faberge, “Apparatus for recording the number of bees leaving and entering a
hive,” J. Sci. Instr., vol. 20, no. 2, pp. 28–311, 1943.
[122] H. G. Spangler, “Photoelectrical counting of outgoing and incoming honey-bees,”
J. Econ. Entomol., vol. 62, no. 5, pp. 1183–1184, 1969.
[123] R. M. Burrill and A. Dietz, “The response of honey bees to variations in solar
radiation and temperature,” Apidologie, vol. 12, no. 4, pp. 319–328, 1981.
[124] C. Liu, J. Leonard, and J. J. Feddes, “Automated monitoring of ﬂight activity ata beehive entrance using infrared light sensors,” J. Apic. Res., vol. 29, no. 1, pp.
20–27, 1990.
[125] R. G. Danka and L. D. Beaman, “Flight activity of USDA-ARS Russian honeybees (Hymenoptera: Apidae) during pollination of lowbush blueberries in Maine,”
J. Econ. Entomol., vol. 100, no. 2, pp. 267–272, 2007.
[126] S. S. Schneider, J. A. Stamps, and N. E. Gary, “The vibration dance of the honeybee. I. Communication regulating foraging on two time scales,” Anim. Behav., vol.
34, no. 2, pp. 377–385, 1986.

118

BIBLIOGRAPHY[127] S. Streit, F. Bock, C. W. W. Pirk, and J. Tautz, “Automatic life-long monitoring ofindividual insect behaviour now possible,” Zoology, vol. 106, no. 3, pp. 169–171,
2003.
[128] J. A. Shaw, N. L. Seldomridge, D. L. Dunkle, P. W. Nugent, L. H. Spangler, J.
J. Bromenshenk, C. B. Henderson, J. H. Churnside, and J. J. Wilson, “Polarization
lidar measurements of honey bees in ﬂight for locating land mines,” Opt. Express.,
vol. 13, no. 15, pp. 5853–5863, 2005.
[129] D. S. Hoﬀman, A. R. Nehrir, K. S. Repasky, J. A. Shaw, and J. L. Carlsten, “Rangeresolved optical detection of honeybees by use of wing-beat modulation of scattered
light for locating land mines,” Appl. Opt., vol. 46, no. 15, pp. 3007–3012, 2007.
[130] E. S. Carlsten, G. R. Wicks, K. S. Repasky, J. K. Carlsten, J. J. Bromenshenk, andC. B. Henderson, “Field demonstration of a scanning lidar and detection algorithm
for spatially mapping honeybees for biological detection of land mines,” Appl. Opt.,
vol. 50, no. 14, pp. 2112-2123, 2011.
[131] V. Estivill-Castro, D. C. Lattin, F. Suraweera, and V. Vithanage, “Tracking bees –a 3D, outdoor small object environment,” Proc. Int. Conf. Image Process., 2003.
[132] T. Landgraf and R. Rojas. “Tracking honey bee dances from sparse optical ﬂow
ﬁelds,” Tech. Rep. B 07-11, Freie Universitat Berlin, 2007.
[133] J. M. Campbell, D. C. Dahn, and D. A. J. Ryan, “Capacitance-based sensor formonitoring bees passing through a tunnel,” Meas. Science Technol., vol. 16, no. 12,
pp. 2503, 2005.
[134] J. Santoyo, W. Azarcoya, M. Valencia, A. Torres, and J. Salas, “Frequency analysis
of a bumblebee (bombus impatiens) wingbeat,” Patt. Anal. Appl., vol. 19, no. 2,
pp.487–493, 2016.

BIBLIOGRAPHY

119[135] A. Decourtye, J. Devillers, P. Aupinel, F. Brun, C. Bagnis, J. Fourrier, and M.
Gauthier, “Honeybee tracking with microchips: a new methodology to measure theeﬀects of pesticides,” Ecotoxicology, vol. 20, no. 2, pp. 429–437, 2011.
[136] J. Campbell, L. Mummert, and R. Sukthanker, “Video monitoring of honey beecolonies at the hive entrance,” Proc. Workshop Visual Observ. Anal. Anim. Insect
Behav., vol. 7, 2008.
[137] T. Kimura, M. Ohashi, K. Crailsheim, T. Schmickl, R. Okada, G. Radspieler, andH. Ikeno, “Development of a new method to track multiple honey bees with complex
behaviors on a ﬂat laboratory arena,” PLoS One, vol. 9, no. 1, pp. 1–12, 2014.
[138] D. J. Kale, R. Tashakkori, and R. M. Parry, “Automated beehive surveillance using
computer vision,” Proc. IEEE SoutheastCon, 2015.
[139] G. J. Tu, M. K. Hansen, P. Kryger, and P. Ahrendt, “Automatic behaviour analysis
system for honeybees using computer vision,” Comput. Electron. Agric., vol. 122,
pp. 10–18, 2016.
[140] G. Chiron, P. Gomez-Kramer, and M. Menard, “Detecting and tracking honeybees
in 3D at the beehive entrance using stereo vision,” EURASIP J. Image Video Process., vol. 1, pp. 1–17, 2013.
[141] E. Caetano, S. Silva, and J. Bateira, “A vision system for vibration monitoring ofcivil engineering structures,” Exp. Techn., vol. 35, no. 4, pp. 74–82, 2011.
[142] H. G. Maas and U. Hampel, “Photogrammetric techniques in civil engineering
material testing and structure monitoring,” Photogram. Eng. Remote Sens., vol. 72,
no. 1, pp. 39–45, 2006.
[143] J. G. Chen, N. Wadhwa, F. Durand, W. T. Freeman, and O. Buyukozturk, “De-120

BIBLIOGRAPHY

velopments with motion magniﬁcation for structural modal identiﬁcation throughcamera video,” Dyn. Civil Struct., vol. 2, pp. 49–57, 2015.
[144] J. Lohscheller, U. Eysholdt, H. Toy, and H. Dollinger, “Phonovibrography: Mapping high-speed movies of vocal fold vibrations into 2-D diagrams for visualizing
and analyzing the underlying laryngeal dynamics,” IEEE Trans. Med. Imaging, vol.
27, no. 3, pp. 300–309, 2008.
[145] D. D. Mehta, D. D. Deliyski, T. F. Quatieri, and R. E. Hillman, “Automatedmeasurement of vocal fold vibratory asymmetry from high-speed videoendoscopy
recordings,” J. Speech Lang. Hear. Res. vol. 54, no. 1, pp. 47–54, 2011.
[146] A. P. Pinheiro, D. E. Stewart, C. D. Maciel, J. C. Pereira, and S. Oliveira, “Analysis of nonlinear dynamics of vocal folds using highspeed video observation andbiomechanical modeling,” Digit. Signal Process., vol. 22, no. 2, pp. 304–313, 2012.
[147] K. Shimasaki, M. Jiang, T. Takaki, I. Ishii, and K. Yamamoto, “HFR-video-based
honeybee activity sensing using pixel-level short-time Fourier transform,” Proc.
IEEE Sensors 2018, pp. 812–815, 2018.
[148] K. S. Repasky, J. A. Shaw, R. Scheppele, C. Melton, J. L. Carsten, and L. H. Spangler, “Optical detection of honeybees by use of wing-beat modulation of scatteredlaser light for locating explosives and land mines,” Appl. Opt., vol. 45, no. 8, pp.
1839–1843, 2006.
[149] A. M. Wenner, “The ﬂight speed of honeybees: a quantitative approach,” J. Apic.
Res., vol. 2, pp. 25–32, 1963.
[150] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi, “You only look once: uniﬁed,
real-time object detection,” IEEE Int. Conf. Comput. Vision and Pattern Recognit.,
pp. 779–788, 2016.

BIBLIOGRAPHY

121[151] G. Levi, and T. Hassner, “Age and gender classiﬁcation using convolutional neuralnetworks,” IEEE Int. Conf. Comput. Vision and Pattern Recognit., pp. 34–42, 2015.

Acknowledgment
Firstly, I would like to thank my advisor, Prof. Idaku Ishii, who gave me an
opportunity to join the robotics engineering laboratory and has guided me through these
six years. He has taught me lot of skills and attitudes for academic research in aspects ofa new research consideration, a research process method, a research summarization suchas presentations and papers. I want to express my gratitude again to Prof. Idaku Ishii.
Additionally, I would like to express my gratitude to Dr. Takeshi Takaki, Dr.
Qingyi Gu, Dr. Taku Senoo, Dr. Mingjun Jiang, Dr. Shaopeng Hu, Dr. Sushil
Raut, and Dr. Michiaki Inoue. They helped me overcome the unfamiliarity with a
new experimental environment when I joined our laboratory. I would like to express myheartfelt gratitude to Ms. Yukari Kaneyuki, and Ms. Michiko Kanzaki (educational
administrator), and Ms. Rumi Horiuchi and Ms. Arisa Tomura (laboratory secretary).
They were my most reliable staﬀ in our institution; I received thoughtful attention bothin my study. I want to express my sincere thanks to the bachelor, master, and doctoral
students in Robotics Laboratory for their help in my life and research. I am gratefulto the teachers, clerical staﬀ, and doctoral students, especially Ms. Diana Khan in the
TAOYAKA Program. They had the opportunity to experience many ﬁeldwork and events
that general doctoral students do not have as well. I had come across valuable experiencethrough onsite course rotation when I was a student in the TAOYAKA program. I will justsubscribe to my point for those experiences regarding onsite projects and training frompoints of the cultural, social, and technical aspects.

123

124

ACKNOWLEDGMENT

Finally, I want to express my profound gratitude to my family for their support andwarm encouragement until now throughout my life.
July, 2020
Kohei Shimasaki