IIEEJ Transactions on
IIEEJ Transactions on
Image Electronics and Visual Computing

Image Electronics andVisual Computing

Vol. 10, No. 2 2022

I.I.E.E.J

Editorial Committee of IIEEJ
Editor in Chief
Osamu UCHIDA (Tokai University)

Vice Editors in Chief
Naoki KOBAYASHI (Saitama Medical University)
Yuriko TAKESHIMA (Tokyo University of Technology)

Advisory Board
Yasuhiko YASUDA (Waseda University Emeritus)
Hideyoshi TOMINAGA (Waseda University Emeritus)
Kazumi KOMIYA (Kanagawa Institute of Technology)
Fumitaka ONO (Tokyo Polytechnic University Emeritus)
Yoshinori HATORI (Tokyo Institute of Technology)
Mitsuji MATSUMOTO (Waseda University Emeritus)
Kiyoshi TANAKA (Shinshu University)
Shigeo KATO (Utsunomiya University Emeritus)
Mei KODAMA (Hiroshima University)

Editors
Yoshinori ARAI (Tokyo Polytechnic University)
Chee Seng CHAN (University of Malaya)
Naiwala P. CHANDRASIRI (Kogakuin University)
Chinthaka PREMACHANDRA (Shibaura Institute of Technology)
Makoto FUJISAWA(University of Tsukuba)
Issei FUJISHIRO (Keio University)
Kazuhiko HAMAMOTO (Tokai University)
Madoka HASEGAWA (Utsunomiya University)
Ryosuke HIGASHIKATA (FUJIFILM Business Innovation Corp.)
Yuki IGARASHI (Ochanomizu University)
Tomokazu ISHIKAWA（Toyo University）
Masahiro ISHIKAWA （Saitama Medical University）
Naoto KAWAMURA (Canon OB)
Shunichi KIMURA (FUJIFILM Business Innovation Corp.)
Shoji KURAKAKE (NTT DOCOMO)
Kazuto KAMIKURA (Tokyo Polytechnic University)
Takashi KANAI (The University of Tokyo)
Tetsuro KUGE (NHK Engineering System, Inc.)
Koji MAKITA (Canon Inc.)
Tomoaki MORIYA (Tokyo Denki University)
Paramesran RAVEENDRAN (University of Malaya)
Kaisei SAKURAI (DWANGO Co., Ltd.)
Koki SATO (Shonan Institute of Technology)
Syuhei SATO（Hosei University）
Masanori SEKINO (FUJIFILM Business Innovation Corp.)
Kazuma SHINODA (Utsunomiya University)
Mikio SHINYA (Toho University)
Shinichi SHIRAKAWA (Aoyama Gakuin University)
Kenichi TANAKA (Nagasaki Institute of Applied Science)
Yukihiro TSUBOSHITA (Fuji Xerox Co., Ltd.)
Daisuke TSUDA (Shinshu University)
Masahiro TOYOURA (University of Yamanashi)
Kazutake UEHIRA (Kanagawa Institute of Technology)
Yuichiro YAMADA (Genesis Commerce Co.,Ltd.)
Norimasa YOSHIDA (Nihon University)
Toshihiko WAKAHARA (Fukuoka Institute of Technology OB)
Kok Sheik WONG (Monash University Malaysia)

Reviewer
Hernan AGUIRRE (Shinshu University)
Kenichi ARAKAWA (NTT Advanced Technology Corporation)
Shoichi ARAKI (Panasonic Corporation)
Tomohiko ARIKAWA (NTT Electronics Corporation)
Yue BAO (Tokyo City University)
Nordin BIN RAMLI (MIMOS Berhad)
Yoong Choon CHANG (Multimedia University)
Robin Bing-Yu CHEN (National Taiwan University)
Kiyonari FUKUE (Tokai University)
Mochamad HARIADI (Sepuluh Nopember Institute of Technology)
Masaki HAYASHI (UPPSALA University)
Takahiro HONGU (NEC Engineering Ltd.)

Yuukou HORITA (University of Toyama)
Takayuki ITO (Ochanomizu University)
Masahiro IWAHASHI (Nagaoka University of Technology)
Munetoshi IWAKIRI (National Defense Academy of Japan)
Yoshihiro KANAMORI (University of Tsukuba)
Shun-ichi KANEKO (Hokkaido University)
Yousun KANG (Tokyo Polytechnic University)
Pizzanu KANONGCHAIYOS (Chulalongkorn University)
Hidetoshi KATSUMA (Tama Art University OB)
Masaki KITAGO (Canon Inc.)
Akiyuki KODATE (Tsuda College)
Hideki KOMAGATA (Saitama Medical University)
Yushi KOMACHI (Kokushikan University)
Toshihiro KOMMA (Tokyo Metropolitan University)
Tsuneya KURIHARA (Hitachi, Ltd.)
Toshiharu KUROSAWA (Matsushita Electric Industrial Co., Ltd. OB)
Kazufumi KANEDA (Hiroshima University)
Itaru KANEKO (Tokyo Polytechnic University)
Teck Chaw LING (University of Malaya)
Chu Kiong LOO (University of Malaya) F
Xiaoyang MAO (University of Yamanashi)
Koichi MATSUDA (Iwate Prefectural University)
Makoto MATSUKI (NTT Quaris Corporation OB)
Takeshi MITA (Toshiba Corporation)
Hideki MITSUMINE (NHK Science & Technology Research Laboratories)
Shigeo MORISHIMA (Waseda University)
Kouichi MUTSUURA (Shinsyu University)
Yasuhiro NAKAMURA (National Defense Academy of Japan)
Kazuhiro NOTOMI (Kanagawa Institute of Technology)
Takao ONOYE (Osaka University)
Hidefumi OSAWA (Canon Inc.)
Keat Keong PHANG (University of Malaya)
Fumihiko SAITO (Gifu University)
Takafumi SAITO (Tokyo University of Agriculture and Technology)
Tsuyoshi SAITO (Tokyo Institute of Technology)
Machiko SATO (Tokyo Polytechnic University Emeritus)
Takayoshi SEMASA (Mitsubishi Electric Corp. OB)
Kaoru SEZAKI (The University of Tokyo)
Jun SHIMAMURA (NTT)Tomoyoshi SHIMOBABA (Chiba University)
Katsuyuki SHINOHARA (Kogakuin University)
Keiichiro SHIRAI (Shinshu University)
Eiji SUGISAKI (N-Design Inc. (Japan), DawnPurple Inc.(Philippines))
Kunihiko TAKANO (Tokyo Metropolitan College of Industrial Technology)
Yoshiki TANAKA (Chukyo Medical Corporation)
Youichi TAKASHIMA (NTT)Tokiichiro TAKAHASHI (Tokyo Denki University)
Yukinobu TANIGUCHI (NTT)Nobuji TETSUTANI (Tokyo Denki University)
Hiroyuki TSUJI (Kanagawa Institute of Technology)
Hiroko YABUSHITA (NTT)
Masahiro YANAGIHARA (KDDI R&D Laboratories)
Ryuji YAMAZAKI (Panasonic Corporation)IIEEJ OfficeOsamu UKIGAYA
Rieko FUKUSHIMA
Kyoko HONDA

Contact Information
The Institute of Image Electronics Engineers of Japan (IIEEJ）
3-35-4-101, Arakawa, Arakawa-ku, Tokyo 116-0002, Japan
Tel : +81-3-5615-2893 Fax : +81-3-5615-2894
E-mail : hensyu@iieej.org
http://www.iieej.org/ (in Japanese)
http://www.iieej.org/en/ (in English)
http://www.facebook.com/IIEEJ (in Japanese)
http://www.facebook.com/IIEEJ.E (in English)

IIEEJ Transactions on
Image Electronics and Visual Computing
Vol.10 No.2 December 2022
CONTENTS

Contributed Papers
141

Audio-Visual Judgement of In-Play or Out-of-Play in Volleyball
Match Videos

Reina ISHIKAWA, Ryo HACHIUMA, Hideo SAITO,
Ryo FUJIWARA, Shoji YACHIDA

150

Multi-Sized Particle Sampling Method Based on Porosity
Optimization in 2D Space

Xu WANG, Makoto FUJISAWA, Masahiko MIKAWA

162

Real-Time Rendering of Oil Film with Flexible Properties

Subroto Prasetyo HUDIONO, Tomoya ITO,
Yuriko TAKESHIMA, Tsukasa KIKUCHI

170

Eyeliner Looked “Attractive” by Eyelids Shape Using Eyeliner
ModelsRena OKURI, Shuhei KODAMA, Tokiichiro TAKAHASHI

Guide for Authors
181

Guidance for Paper Submission

Published two times a year by the Institute of Image Electronics Engineers of Japan （IIEEJ）
3-35-4-101，Arakawa，Arakawa-ku，Tokyo 116-0002，Japan
Tel：+81-5615-2893 Fax：+81-5615-2894 E-mail：hensyu@iieej.org http://www.iieej.org/

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）
IIEEJ

Paper -- Special Issue on Journal Track Papers in IEVC2021 --Audio-Visual Judgement of In-Play or Out-of-Play in Volleyball Match Videos
Reina ISHIKAWA† , Ryo HACHIUMA† , Hideo SAITO† (Member), Ryo FUJIWARA†† , Shoji YACHIDA††† Keio University, Japan , †† NEC, Japan<Summary> Utilizing deep learning in sports match analysis is very promising for eﬃciency, and there
have been many studies focusing on volleyball. However, previous studies have used only the frames where
the ball is alive (i.e., in-play frames), and recent data analysis methodologies rely on visual data, which arehighly dependent on the position and angle of the camera relative to the court. As a result, these methods
require a large dataset of images taken from various angles to improve accuracy. To eﬃciently extract andanalyze the play as it is happening, we propose a model that distinguishes between in-play and out-of-playby combining visual data and audio data of volleyball match videos using late fusion. To investigate theeﬀectiveness of the proposed model, we perform Grad-CAM visualization to determine which pixels the
proposed model is focusing on.
Keywords: volleyball analysis, binary classiﬁcation, audio-visual1.

Introduction

Sports analysis is beneﬁcial to understand thestrengths and conditions of a team and each player
therein. By analyzing and quantifying the players’ conditions in real-time during a match, the coach can give
appropriate advice to the players or replace them withplayers who are better suited to the situation. Furthermore, by reviewing the analysis after a match, players
can identify their weaknesses and use that information
to develop eﬀective training plans and new team tactics.
Conventionally, sports analysis had been implemented
manually by specialists who are referred to as sports analysts. However, neural networks have appeared, andthe idea of integrating sports analysis with neural networks has the potential power to remove conventional
manual eﬀort and limitations. This integration also enables a more eﬃcient and detailed analysis. In particular, the automation of analysis in team sports suchas volleyball and basketball, which are diﬃcult to analyze manually, is in great demand given its ability to
improve the eﬃciency and accuracy of the work. For
example, individual or group action recognition 1)–4)andplayer detection and tracking studies 5),6)provide an es-

sential foundation for analysis. Azar et al.1)proposed amethod that undertakes both individual and group action recognition tasks to reﬂect individual-level actions
in group-action recognition.
The problem with the conventional studies referenced
above is that they use datasets that are supposed tocontain in play data, so it is assumed that the out-ofplay frames are removed beforehand. In this paper, byin play, we mean that the ball is alive, that is, it isin the time period when a referee blows the whistle tostart the play to when the next whistle is blown. On
the other hand, out of play means that the ball is dead,which refers to a period other than in-play. As volleyballis characterized by the fact that the out-of-play time
is longer than the in-play time throughout the match,eliminating out-of-play time can lead to more eﬃcient
analysis. In this context, the automatic determinationof whether the ball is in play is highly important forthe downstream volleyball video analysis tasks, such as
action recognition and ball/player tracking.
Therefore, we present a method that categorizes thescene as in play or out of play from the volleyball recordings. As our method eliminates redundant moments

141

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）…

Image
Subnetwork…

In play

FusionNetwork

Grayscale framesOR
Out of play

Audio
Subnetwork

Original film

Cochleogram

Fig. 1 Overview of the proposed method(out-of-play frames) from the recorded video, it has thepotential to speed up the volleyball analytics. By eliminating out-of-play scenes that are not useful in analytics,
we expect that team staﬀ, including players and coaches,
can review their rallies in a shorter time. This elimination can let them smoothly analyze their states, such as
their strengths or weakness. In addition, when we apply
another process to the match video for further analysis,for example, using machine learning, eliminating redundant periods can reduce processing time.
We note that we are targeting volleyball staﬀ, including volleyball analysts who are specialists in volleyballmatch analysis and coaches. Therefore in this paper, we
assume to use raw videos taken by team members, andwe do not use broadcast videos with replay scenes oroverlaid graphics, such as scores and subtitles.
Moreover, most conventional prediction models1)–4)
make decisions based only on visual data, and existing datasets contain only images in which the camera
angle does not change signiﬁcantly and the sight line isalmost parallel to the centerline. However, such conventional methods cannot cope with signiﬁcant changes
in the viewpoint. Therefore, using only visual data requires a large amount of image data from various angles.
In this study, we propose a method that uses audio datain addition to visual data. The audio data is considered
to be less sensitive to changes in the angle of the viewpoint than the visual data and is thus expected to provide more robust judgment results. As far as we know,
prior to this study, audio data has never been used involleyball analysis.
Our contributions in this paper are as follows:• We propose a framework that uses both audio andvisual data to determine whether a given second in
a volleyball game is in play or not using a convolutional neural network (CNN).
• We conduct experiments comparing the proposed
method to the conventional method using either image or audio data. These experiments demonstrate
how robust the proposed model is against changes
in the seat occupancy rate or the camera position.
• We visualize the image subnetwork of the proposedmodel using Grad-CAM.

2.

Related Work

Humans understand their surrounding environments
in a multi-modal manner, using ﬁve senses: sight, touch,taste, smell, and hearing. When conﬁdence in one
modality is low, we can use the remaining modalities
to make decisions; in other words, they can compensate
for one other. These biological structures have recently
inspired the investigation of multi-modal sensor usage invarious ﬁelds. Here we introduce some previous works
that have used multi-modalities and have focused on
methods that incorporate audio features and visual features, which are used in this paper.
Understanding the ground conditions on which anautonomous driving robot will operate is crucial for
ensuring smooth movement7),8). For example, Zürn
et al.9)proposed a self-supervised semantic segmentation
method that takes advantage of independent audio andvisual feature extraction.
The audio-visual feature is also utilized in the ﬁeld ofsports analysis. For instance, Xiong et al.10) proposed a142

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

highlight extraction method that associates visual object (e.g., baseball catcher and soccer goal post) detection and audio markers. According to the authors of10),
each of the visual markers and audio markers can compensate for the other’s false detection, which is a similar idea to the one we explore here. Khan11) proposed
an audio-visual highlight generalization method focusing on broadcast sports videos. They take advantage
of audio features, such as the cheering of the audience,
to identify exciting points within the video. Similarly,
Grade et al.12) also proposed a sports-type classiﬁcation
method that combines visual motion features and audio feature processed with mel frequency cepstral coeﬃcients (MFCCs) and principal component analysis(PCA). They successfully classify three sports types:
basketball, soccer, and volleyball. The audio-visual usage is similar to ours, but the diﬀerence is the targetof the classiﬁcation. The previous study classiﬁes the
sports type and does not further split a single video.
On the other hand, ours classify one video into many
parts.
Itazuri et al.13)proposed a method to extract rally
scenes by scoring camera transition and sound information. They assume to take the court in a speciﬁc
direction because they are thinking of using the videos
which are open to the public, and the users of the gameare the general public. In contrast, our method targets
sports analysts in a team who use videos taken by teammembers from their bench or audience seats. One of ourcontributions is robustness to the change of the camera
position due to the usage of both visual information andaudio information. Therefore, our study is uniquely positioned in terms of utilizing audio-visual data to recognize whether segments are in play or not in play involleyball matches.

3.

Methodology

In this study, segments are deﬁned as one second, andthe goal is to categorize each segment as either in play orout of play (see Fig. 1). For visual data, 10 frames are
sampled uniformly from the target segment, and eachsegment is resized and converted into 144 × 144 × 10
grayscale images. Then, the visual input is created byrandomly cropping a 128 × 128 image from the original grayscale frame and combining the 10 grayscale images. To process a one-dimensional audio signal into a
two-dimensional input signal, we adopted cochleogram.
Though there are other audio pre-processing methods,
such as Mel-frequency cepstral coeﬃcients (MFCCs),
we suppose pitch information kept with cochleogram isessential to diﬀerentiate the sound of hitting the balland the sounds of cheers, for example, applause. For
this purpose, each audio segment is converted into twodimensional data of 64 × 64 using cochleogram and used
as input.
The image and audio subnetworks are constructedbased on the VGG1615). However, the number of layers are changed from the VGG16 according to the input
size, and the fully connected layers are removed. Eachkernel size is 3 × 3, and a max-pooling layer follows a setof convolutional layers. The image subnetwork consistsof six layers, and the audio subnetwork consists of ﬁve
layers. In the end, both subnetworks output a feature
on a 256 × 2 × 2 scale. Each feature output from theimage and audio subnetworks is reshaped into a 1024dimensional vector, and those are added and then input
into the fusion network, which consists of 3 fully connected layers with a 0.5 dropout rate and a ReLU activation function. The output of the fusion network is
a two-dimensional vector. Each element of the vector
indicates the score for each class (i.e., in play or out ofplay), and the ﬁnal judgment is made based on whichclass has the higher score value.

4.
4.1

Experiment

DatasetIn this study, we created a dataset from eight variety volleyball match videos to evaluate the proposed
method. A sample image from each video is shown in
Fig. 2. All of the volleyball match videos used in thispaper were captured by a staﬀ of the volleyball teamplaying from a static position. Among the videos takenfrom angles in #1–#8 of Fig. 2, a video from a particularly diﬀerent angle is displayed in #4. The seatoccupancy rates were 100% in the 1st , 2nd , 3rd , and 4thvideos, 50% in the 5th and 6th videos, and almost 0%

143

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

#1

#2

#3

#4

#5

#6

#7

#8

Fig. 2 Sample images from each video
except for people aﬃliated with the team in the 7th and8th videos. We manually annotate the approximatedoccupancy rates for videos.
Each video is segmented into 1-second intervals, andthe image and audio for each segment are preprocessed
as explained in Sec. 3. before they are put into each subnetwork. Table 1 shows both the length and the number of Segments in each video. To create the groundtruth labels, we manually labeled only the 4th video in
Fig. 1. As for the remaining videos, we utilized thelabels created by volleyball analysts in real-time while
the match was going on. The part we want to extract as
in-play is the period from when the 1st referee blows awhistle of starting a rally to when the referee blows thenext whistle. However, in volleyball, there is a rule thata server has to carry out a service hit within 8 seconds
after the 1st referee’s whistle. In other words, there may
be an up to 8 seconds gap between the starting whistle
and the time when the server hits the ball. Because analysts record the approximate time when a player tosses
or releases the ball and hits it, their recorded time may
diﬀer from the actual time of the whistle. To ﬁll in thegaps as much as possible, we manually labeled the threevideos based on the refrees’ whistles, which we assumed
was the correct labeling, and we decided to take a margin to ensure the highest level of the accuracy of thelabeling based on the analysts’ data. Adding margins of
2 seconds for the start of in play and +5 seconds for theend resulted in a 95.09% match between manual labelingand analysts’ labeling. The total number of in-play seg-

Table 1 The seat occupancy rate and the length ofand number of segments in each video
Video
1
2
3
4
5
6
7
8

Seat occupancy
rate (%)
100
100
100
100
50
50
0
0

Length# Segments

01:56:31
01:41:23
02:30:52
00:26:39
02:25:53
01:48:37
02:18:49
02:07:13

6991
6083
9052
1599
8753
6517
8329
7633

ments is 19,996, and the number of out-of-play segments
is 34,961.
4.2

Dataset splitting

We divided the dataset we generated into separate
training, validation, and test datasets with two versions.
We prepared a heavy version and a light version for eachsplitting method according to the seat occupancy rate
shown in Table 1. Because the heavy version was utilized
to examine the impact of the seat occupancy rate, it
includes all eight recordings. In contrast, since the lightversion was utilized to examine the impact of diﬀerencesin the court orientation, it includes the recordings thathave a 100% seat occupancy rate.
Random split We randomly assigned 80% of thedataset to be part of a training dataset, 10% to be partof a validation dataset, and 10% to be part of a test
dataset.
Time series split A split for training is made suchthat 80% of the time series video is allocated from the144

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

beginning. The test dataset is from the next 10%, andthe validation dataset is from the remaining 10%. The
number of data in each dataset is almost the same as in
the random split.
Sequence split We divided the dataset by video unit
as shown in Table 1. First, we separated all the datasets
into three groups according to the seat occupancy rate
shown in the same table: 100%, 50%, and 0%. Then,each of the heavy and light versions was divided in thefollowing ways:• To investigate the impact of the seat occupancyrate, we implement leave-one-out classiﬁcations,that is, one of the groups is assigned to the test
dataset, and the remaining groups are assigned tothe training and validation datasets.
• The 1st , 2nd , and 3rd videos, which have similar
court orientations, were assigned to the training
and validation datasets. The data from the 4thvideo, which has the most variation in court orientation, was assigned to the test dataset. Thisdata splitting is used to evaluate the eﬀectivenessof audio data when the camera position relative tothe court in the test video was diﬀerent from thetraining dataset (see Subsec. 4.5.2).
Random splitting is the easiest split method for themodel to learn the data distribution in each video, andit works as a benchmark. When the accuracy is high using the random split method and drops when using othersplitting methods, the model is overﬁtting to short-timefeatures. In this case, the model cannot classify the input data, even if it has been trained using the video withthe same court orientation, and the method is neverpractically applicable. When the accuracy is high inthe time series split and drops in the sequence split, themethod can be applied only to other match videos using
a similar court orientation. Our ﬁnal goal is to achieve
high accuracy even for unknown match videos, that is,
to improve the accuracy by sequence segmentation.
4.3

Training details

Because the number of data per class is unbalanced,
we adopt weighted cross-entropy for the loss function
whose weights are hyperparameters. We set the weightof the in-play class to 3 and the weight of the out-ofplay class to 1. We train our model using the Adam
optimizer16) with a learning rate of 1e-6. The iteration
number in the training time is set to 50 epochs for thelight dataset (i.e., the dataset that uses the 1st to the
4th videos) and 60 epochs for heavy dataset (i.e., thedataset that uses the 1st to the 8th videos).
We used the model with the highest F-value against
the validation data for testing. All of our implementations were based on the Pytorch framework,14) and all
the experiments were performed on a computer with anIntel(R) Xeon(R) W-3235 CPU and a Titan RTX GPU,
64GB memory.
4.4

Evaluation metricWe evaluated our model by calculating the accuracy,
precision, recall, and F-score for each of the splitting
methods, with positive indicating that the ball was alive
and negative meaning that the ball was dead. The precision was determined based on the percentage of segments with positive labels among the segments the inference network judged as positive. The recall was determined by the percentage of the segments that theinference network judged as positive among those withpositive labels. The F-score is the harmonic mean of
the precision and the recall rates, which are deﬁned asfollows:
2

F =

.
1
1
recall + precision

4.5

Evaluation(1)

To evaluate the eﬀectiveness of our method, we created the visual model that uses only visual data and theaudio model that uses only audio data, and we compared
our method with those two models, respectively. The visual model receives only visual input and is similar to
the image subnetwork. It has fully connected layers forclassifying into in play or out of play in addition to theimage subnetwork structure. Similarly, the audio model
receives only audio input and has fully connected layers
for classiﬁcation after the audio subnetwork structure.
4.5.1

Seat occupancy rate

As the seat occupancy rate increases, the cheers fromthe audience get louder, which the audience cheering can

145

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Table 2 Comparison between the audio model and our model using datasets with diﬀerent seat occupancy rate
Dataset Split
RandomTime Series
Sequence(test with 100% seat occupancy rate videos)
Sequence(test with 50% seat occupancy rate videos)
Sequence(test with 0% seat occupancy rate videos)

Model
Audio
Ours
Audio
Ours
Audio
Ours
Audio
Ours
Audio
Ours

aﬀect the result of the classiﬁcation task. Therefore, toinvestigate robustness against changes in the seat occupancy rate, we divided the dataset into three groups according to the seat occupancy rates, 100%, 50%, and 0%,
and implemented leave-one-out training among thesethree groups. This is referred to as a heavy sequence
split. This data split method requires the model to infer over the data with unknown seat occupancy rates.
For instance, when we train our model using the data of
100% and 50% seat occupancy rates, we test our modelusing the data with a 0% seat occupancy rate. We also
compared with those acquired using the heavy random
split dataset and the heavy time series split dataset.
The results are shown in Table 2. With random split
or time-series data split, our model surpassed the audio
model in terms of precision, recall, and F-score. With
sequence split, when the seat occupancy rate was 100%,
the audio model outperforms ours in terms of both accuracy and F-score. This seems to result from the audience cheering signiﬁcantly impacting the in-play andout-of-play classiﬁcations, and under the full occupancy
rate condition, using only audio data can yield high accuracy. However, as the seat occupancy rate decreases,
the precision rate of the audio model rapidly decreases
because the audience cheering gets low. On the whole,
our model demonstrates its robustness to the decrease of
the seat occupancy rate. This result reveals that using
an image subnetwork is eﬀective.
4.5.2

Court orientation

In this subsection, we explore the eﬀect of court orientation on our model using the 100% seat occupancy
dataset, that is, the light version dataset. Among the

Accuracy↑
0.804
0.822
0.792
0.748
0.862
0.807
0.789
0.856
0.663
0.719

Precision↑
0.685
0.695
0.668
0.741
0.758
0.736
0.662
0.794
0.524
0.576

Recall↑
0.845
0.903
0.850
0.870
0.910
0.732
0.854
0.811
0.901
0.889

F-score↑
0.757
0.785
0.748
0.800
0.827
0.734
0.746
0.803
0.663
0.699

four videos, the 4th one has the most variation in court
orientation of all the videos. We can observe from theresult of Table 3 that with the sequence split dataset,
where there were changes to the court orientation of thetest data as compared to the random split or time series split datasets, the F-score and the accuracy dropped
when we used the visual model. On the other hand, theproposed method, which uses both audio data and image data, resulted in better F-scores and accuracy withthe sequence split than the visual model. Therefore, we
can say that the proposed method is more robust tochanges in the camera view than the visual model, andthis is owing to the usage of the audio subnetwork.
4.5.3

Visualization with GRAD-CAM

Subsection 4.5.2 demonstrated that using only visual
data on its own may lead to erroneous decisions when
the court orientation changes signiﬁcantly, and the visual data may provide no information about the court
orientation at all. With this in mind, in this subsection, we examine whether the visual network can capture information about the court and help the model
make decisions. Speciﬁcally, we used Grad-CAM17) tovisualize the model trained on the sequence split dataset
and examine which pixels are of interest in the visual input. Grad-CAM observes the gradients of inputs in each
CNN layer to visualize which part of the input image issigniﬁcant for predictions.
We show the visualization results in Fig. 3. From left
to right, each image corresponds to the 1st to 5th layer inthe proposed image subnetwork. Figure 3(a) represents
the visualization of the layers by weighting the activation map based on the gradient of the in-play score for146

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Table 3 Evaluation using diﬀerent court orientation video data
Dataset split
RandomTime Series
Sequence

Model
VisualOurs
VisualOurs
VisualOurs

Accuracy↑
0.705
0.796
0.739
0.811
0.558
0.736

(a) In-play visualization based on the in-play score

Precision↑
0.543
0.640
0.569
0.657
0.507
0.704

Recall↑
0.934
0.937
0.943
0.927
0.729
0.716

F-score↑
0.686
0.761
0.710
0.769
0.598
0.710

(b) Out-of-play visualization based on the out-of-play score

Fig. 3 Visualization samples created with Grad-CAM
an in-play image, and Fig. 3(b) shows the visualizationof the layers by weighting the activation map based on
the gradient of the out-of-play score for an out-of-play
image. For each image group in (a) and (b), the upper result was generated using an image sample fromthe training data, and the lower one was generated using a sample from the test data. The ﬁgure shows thatthe inference network focuses on the players inside thecourt when drawing inferences from in-play images; In
contrast, the network focuses on the unoccupied partof the court or the spectators when drawing inferencesfrom out-of-play ones.
These results derive some considerations; To be detected as in-play, player positions can be informative;
To be judged as out-of-play, the unoccupied part of thecourt or the move of the audience outside the court canbe signiﬁcant. Though the camera positions diﬀer between the training dataset and the test dataset, we cansee the network can detect the court position and theplayers inside the court. Because the featured part ofthe input image is the same regardless of the court orientation, the inference network can work, in other words,
distinguish in-play or out-of-play without regard to thecourt orientation. Furthermore, since the attention areas are distinctly diﬀerent for in-play versus out-of-play,
image data is considered informative for in-play judg-

ment.

5.

Discussion

Though our model demonstrated better results than
audio model when the seat occupancy rate decreases,
our model scored the lowest accuracy when using 0%
seat occupancy rate videos. We assume this is because
the cheering by the audience has a powerful impact on
the in-play and the out-of-play classiﬁcation. However,
the most informative and invariable sounds are the referees’ whistles and the sound of the plays inside thecourt for conducting the classiﬁcation. Therefore, reducing the impact of the audience voices or extracting
the whistles of the referees is expected to improve therobustness of the decrease in seat occupancy rate. Exploring these eﬀects can be our future work.
In addition, the main limitation of this method is thelabeling system. The labels for training are created byvolleyball analysts or by team members during the volleyball match or created by the authors based on thewhistle of the chief referee. The ground-truth boundaries between in-play and out-of-play created by analysts depend on the skill and judgment of the analysts.
Moreover, as for the authors labeling, when the soundof the whistle is blocked by cheers, picking out the whistle sound is diﬃcult, so we need to label referencing the147IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

move of the referees. For the above reasons, in bothlabeling methods, the boundaries between in-play scene
and out-of-play scene can be blurred, and sometimes it
is hard to deﬁne absolute boundaries. Therefore, we
suppose the ground-truth labels can have a slight error,
and it is diﬃcult to achieve 100% accuracy.

6.

ers, Materials & Continua, Vol.58, No.3, pp. 625–639 (2019).
7) A. Kurobe, Y. Nakajima, K. Kitani, H. Saito: “Audio-VisualSelf-Supervised Terrain Type Recognition for Ground Mobile
Platforms,” Journal of IEEE Access, Vol. 9, pp. 29970-29979(2021).
8) R. Ishikawa, R. Hachiuma, H. Saito:“Self-Supervised AudioVisual Feature Learning for Single-Modal Incremental Terrain Type Clustering,” Journal of IEEE Access, Vol. 9, pp.
64346-64357 (2021).

Conclusion

In this paper, we presented an audio-visual binaryclassiﬁcation method for determining whether a given
second in a volleyball match video is in or out of play.
Though the proposed method, which utilizes both audio
and visual data, is quite simple, we showed it is robust to
diﬀerences in the environment for each video since audio
data and visual data can compensate for each other. Because we are targeting volleyball staﬀ, including analysts
and coaches, and we assume they use the practice or ofﬁcial match videos recorded by themselves, robustness
against the changes in the viewing angle and changes inthe seat occupancy rate is very beneﬁcial.
We used Grad-CAM to verify which parts of the volleyball frame the model prioritizes. The Grad-CAM
visualization revealed that the image subnetwork canrecognize the inside and outside of the court, which is
beneﬁcial for making in-play or out-of-play decisions.

References
1) S. M. Azar,
M. G. Atigh,
A. Nickabadi,
A.
Alahi:“Convolutional Relational Machine for Group Activity Recognition,” Proc. of the IEEE/CVF Conference on
Computer Vision and Pattern Recognition, pp. 7884–7893(2019).
2) J. Wu, L. Wang, L. Wang, J. Guo, G. Wu:“Learning Actor
Relation Graphs for Group Activity Recognition,” Proc. of
the Computing Research Repository, pp. 9956–9966 (2019).
3) M. S. Ibrahim, S. Muralidharan, Z. Deng, A. Vahdat, G.
Mori:“A Hierarchical Deep Temporal Model for Group Activity Recognition,” Proc. of IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 1971–1980 (2016).

9) J. Zürn, W. Burgard, A. Valada:“Self-Supervised Visual
Terrain Classiﬁcation from Unsupervised Acoustic Feature
Learning, ” Journal of IEEE Trans. on Robotics, Vol. 37,
No. 2, pp. 466-481(2021).
10) Z. Xiong, R. Radhakrishnan, A. Divakaran, T. S.
Huang:“Highlights extraction from sports video based on an
audio-visual marker detection framework,” Proc. of the IEEE
International Conference on Multimedia and Expo (2005).
11) A. A. Khan, J. Shao, W. Ali, S. Tumrani:“ContentAware Summarization of Broadcast Sports Videos: An Audio–Visual Feature Extraction Approach,” Journal of Neural
Process Letters, Vol.52, No.3, pp1945–1968 (2020).
12) R. Grade, M. Abou-Zleikha, M. G. Christensen, T. B.
Moeslund:“Audio-Visual Classiﬁcation of Sports Types,”
Proc. of the IEEE International Conference on Computer Vision Workshop, pp. 768-773 (2015).
13) T. Itazuri, T. Fukusato, S. Yamaguchi, S. Morishima:“CourtBased Volleyball Video Summarization Focusing on Rally
Scene,” Proc. of the IEEE Conference on Computer Vision
and Pattern Recognition Workshops, pp. 179-186 (2017).
14) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G.
Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A.
Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejan, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala:“PyTorch: An Imperative Style, High-Performance Deep
Learning Library, ”Proc of the Computing Research Repository (2019).
15) K. Simonyan, A. Zisserman:“Very Deep Convolutional Networks for Large-Scale Image Recognition,” Proc. of the International Conference on Learning Representations (2015).
16) D. P. Kingma, J. Ba: “Adam: A Method for Stochastic Optimization,” Proc. of the International Conference on Learning
Representations (2015).
17) R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D.
Parikh, D. Batra:“Grad-CAM: Visual Explanations from
Deep Networks via Gradient-Based Localization,” Proc. of
the IEEE International Conference on Computer Vision, pp.
618–626 (2017).

(Received August 19, 2022)(Revised October 10, 2022)

4) M. Qi, Y. Wang, J. Qin, A. Li, J. Luo, L. Van Gool:“stagNet: An Attentive Semantic RNN for Group Activity andIn-dividual Action Recognition,” Journal of IEEE Trans. on
Circuitsand Systems for Video Technology, Vol.30, No.2, pp.
549–565 (2020).
5) A. A. Sangüesa, G. Haro, C. Ballester:“Multi-Person Tracking by Multi-Scale Detection in Basketball Scenarios,” Proc.
of the Computing Research Repository (2019).
6) P. K. Santhosh, B. Kaarthick:“An Automated Player Detection and Tracking in Basketball Game,” Journal of Comput-

148

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Reina ISHIKAWA
She received her B.E. degree in informationand computer science from Keio University,
Japan, in 2021. She is currently pursuing
her M.Sc.Eng. degree in science and technology at Keio University.

Ryo

HACHIUMA

He received his B.E. and M.Sc.Eng. degree
in information and computer science from
Keio University, Japan, in 2016 and 2017,
respectively. He received Ph.D. degree inscience and technology with Keio University,
Japan in 2021. He is currently working at
Konicaminolta Japan, as a computer visionresearcher.

Hideo SAITO (Member)He received the Ph.D. degree in electricalengineering from Keio University, Japan, in
1992. Since 1992, he has been on the Facultyof Science and Technology, Keio University.
From 1997 to 1999, he joined the Virtualized Reality Project as a Visiting Researcher
at the Robotics Institute, Carnegie Mellon
University. Since 2006, he has been a Full
Professor at the Department of Information
and Computer Science, Keio University.
Ryo

FUJIWARAHe has been an analyst in the NEC Red
Rockets Volleyball team since 2018, where
he has contributed to the team training.
He served as a support analyst for the
Women’s volleyball team of Japan in theTokyo Olympics 2021. Since 2022, he has
joined the Sports business promotional department of NEC Corporation, where he explores the sports DX area using new technologies and data analysis.

Shoji YACHIDA
He joined NEC Home Electronics, Ltd. in
1988. He joined NEC Corporation in 2000
and is now a senior researcher at the VisualIntelligence Research Laboratories of NEC.
He has engaged the research for the highspeed vision sensing technology and its applications.

149

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）
IIEEJ

Paper

-- Special Issue on Visual Computing --Multi-Sized Particle Sampling Method Based on Porosity
Optimization in 2D Space
Xu WANG† , Makoto FUJISAWA† (Member) , Masahiko MIKAWA†† University of Tsukuba<Summary> This paper introduces a multi-sized particle sampling method within an arbitrary 2D shapeusing power tessellation. Our method aims to improve packing density as to sample as many particles as
possible in a limited area. We propose a porosity-driven optimization technique to ensure no overlap between
particles while increasing the packing density. With such properties, our method is applicable to physicallybased simulations, such as the Discrete Element Method and its related framework. Additionally, thistechnology allows users to set the target particle size distribution by a predesigned cumulative distributionfunction and restrict the errors between 10% and 20% after the optimization. We demonstrate that ourmulti-sized particle sampling algorithms signiﬁcantly improve packing density compared to Poisson disk
sampling and SPH-based blue noise sampling.
Keywords: multi-sized particle sampling, porosity optimization, power tessellation, discrete element method

and 20% after the optimization.

1.

Introduction

Uniform particle sampling is a highly applicable technique in many graphics-related ﬁelds, such as image synthesis, rendering, geometry processing, and physicallybased simulations. Notably, it has played an essential
role in particle-based physics simulations in recent years.
Meanwhile, several adaptive particle (multi-sized) sampling methods1)–3)are also implemented into the ﬂuid simulation. However, it is not well suited for simulating granular materials, such as sand, because slight overlap candestabilize the simulation. Therefore, we aim to achieve
multi-sized particle sampling in arbitrary 2D polygons
and have it available for using in stable physics-basedsand simulations.
This paper proposes a multi-sized particle samplingalgorithm using weighted centroidal power tessellation.
It computes the particle distribution by minimizing theporosity in a pre-deﬁned 2D polygon. Our contributions
are as follows:• We present a novel multi-sized particle sampling
method based on minimum porosity. Our method
ensures no overlapping occurs between particles andeﬀectively exploits the limited space to improve sampling density.
• Our method supports users to give a target particle
size distribution using a predesigned cumulative distribution function and limits the error between 10%150• Multi-sized particles sampled by our method areavailable to the Discrete Element Method (DEM)
and can realize stable physical simulations, whereas
the particles sampled in the comparison experiments
fail to perform stable physical simulations.

2.

Related Work

2.1 Uniform particle samplingBlue noise sampling is a well-known method used incomputer graphics because of its ability to generate random uniform distributions. One of its patterns, Poissondisk sampling, has many rendering and geometry processing applications. Bridson4)proposed a fast Poissondisk sampling method that is eﬀective and more easily
implemented in arbitrary dimensions. While this method
has a wide application, it is diﬃcult to control the sampling density and ensure no overlapping between particles. Traditionally, uniform particle sampling within a 2D
polygon is mainly implemented by dividing the plane as
the uniform grids (Voxelization in 3D5)) and generatinga particle in each uniform grid. However, it will resultin a particle distribution that does not ﬁt the boundary shape. As a solution to this problem, Schechteret al.6)proposed a relaxation-based sampling method tosmoothen the particles on the boundary. Later, Jiang etal.7)applied the Smoothed Particle Hydrodynamics (SPH)
method to achieve blue noise sampling, which can also re-IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 1 Schematic diagram of porosity

alize uniform particle sampling by relaxing the boundary
particles. Although these methods work well for uniform
particle sampling, what we expect as multi-sized particle
sampling and stable sand simulation cannot be fulﬁlled
by simply extending their method.
2.2

Adaptive particle sampling

In order to achieve high-performance SPH-based ﬂuid
simulation, Adams et al.1)presented a dynamic sampling
method, which can dynamically adjust the size of thesampled particles by computing the distance between the
medial axis and ﬂuid surface. Jiang et al.7)integrated
the adaptive particle sampling method and proposed an
SPH-based blue noise sampling to treat multi-sized particles. Winchenbach et al.2)proposed an adaptive incompressible SPH-based ﬂuid simulation method and an improved framework3)adapted to higher volume ratios for
simulating low-viscosity turbulent ﬂows and achieving thehigh-resolution rendering of ﬂuid surfaces. However, sincethese algorithms are all SPH-based, there will be overlapping between the multi-sized particles, which make themulti-sized sand simulation unstable.
2.3

Voronoi diagram and centroidal voronoi
tessellation

The Voronoi diagram has been widely used in thearea of computer graphics. It has many applications inphysics-based simulations of computer animation, suchas power particles8), adaptive staggered power particles9),
power diagram-based high-resolution adaptive liquids10).
Over time, several GPU-based algorithms diagram havebeen developed aiming to speed up the computation of
the 3D Voronoi11)–14). Among these methods, centroidal
Voronoi tessellation (CVT) for reaching a speciﬁc density distribution is also one of the critical topics related
to the Voronoi diagram. CVT can be regarded as oneof the relaxation-based sampling methods. It has twostages of computation: 1. Randomly generates the point
set. 2. Using Lloyd iterations to optimize the loca-

151

tion of the points until convergence. However, directly
using CVT will not cause the sampled points to havebetter properties of blue noise. Balzer et al.15),16)proposed a capacity-constrained Voronoi tessellation, which
can generate high-quality blue noise properties and perfectly ﬁt a known density function. According to the feature that the energy function of CVT has second-ordersmoothness, Liu et al.17)presented a faster convergence
CVT method. Chen et al.18)proposed a parallel algorithm
for CVT on GPU and veriﬁed the algorithm’s performance, which can obtain good blue noise characteristics
and high-speed computing. All these methods describedabove have diﬃculty in extending to multi-sized particle
sampling. On the other hand, CVT was used to study
the hierarchical data visualization presented by Balzer etal.19). They called this approach “Voronoi Treemaps”.
Later, Nocaj et al.20)introduced a faster, more straightforward, and resolution-independent method for computing Voronoi treemaps. While this approach can generate
multi-sized particles, the user needs to pre-deﬁne a large
amount of data since its computation is based on hierarchical data. Besides, sampled results are singular withthe same data.
2.4

Physics based granular material simulation

Granular material simulation-related studies have along history in the computer graphics community21)–25). Numerous works in recent years still focus on sandwater mixing simulation frameworks, such as SPH-based
multi-phase multi-material ﬂows26),27)and material point
method (MPM) based simulations of a mixture of sandand water28),29). Since the particle-particle approach(DEM) could better capture small-scale features at theparticle level, several studies have also employed thismethod to make sandy-like animation, e.g., hybrid grains
30)
and seepage ﬂow framework31).

3.

Porosity Optimization Based ApproachThis section will introduce our porosity optimizationbased method of how to sample multi-sized particles.
3.1

Controllable porosity

Generally, a rock or soil’s porosity is measured by thefraction of the total soil volume occupied by pore space
per unit volume. As shown in Fig.1, the porosity ϕ isdeﬁned as:

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 2 Result of standard gridding approachbased uniform particle sampling
ϕ=1−

Vsoil
Vtotal(1)where Vsoil is the volume of soil particles, and Vtotal
presents a unit soil volume. The high ratio of soil particles in the unit volume will produce minor porosity.
Hence we can improve the packing density within the unit
volume by reducing the porosity and ﬁlling the limited
space with more particles. Thanks to the controllabilityof the porosity, we can obtain higher packing densities
and multi-sized particles by using a power diagram andporosity optimization.
3.2

Power diagram

There are n points S = {s1 , s2 , . . . , sn } (sicalled Voronoi site) in the plane, Voronoi diagram
is made by dividing the plane into n regions R ={R (s1 ) , R (s2 ) , . . . , R (sn )} (each site si ∈ S and its associated region R (si ) are called Voronoi cell) and eachregion has only one point. Given an arbitrary 2D polygon as boundary Ω, Voronoi cell R (si ) can be deﬁned
as:
R (si ) = {p ∈ Ω : ∥p − si ∥ < ∥p − s∥for each s ∈ S − si }(2)

When assigning a weight wi to each Voronoi site si ,we call this Voronoi diagram as Power diagram which isdeﬁned as:
R (si , wi ) = {p ∈ Ω : ∥p − si ∥2 − wi < ∥p − s∥2 − w
for each s ∈ S − si } (3)When the weights wi of all Voronoi site si are equal
to zero, the power diagram is equivalent to a Voronoi
diagram.
In contrast to the Voronoi diagram, the power diagram
allows us to control the shape and size of each Voronoi

152

Fig. 3 A ideal case of the geometric relationshipbetween the Voronoi cells and its weights

cell in a relatively simple way. It means that we canmore readily manipulate the size of the sampled particles and the porosity of the entire space. For controllingthe size of Voronoi cell, there are many studies related tothe capacity-constrained method. In these studies, somemethods combined power diagrams with a false position
method to obtain the diﬀerent sizes of Voronoi cells16)or
integrated with Newton’s method32), L-BFGS method33). However, all these methods have a high computational
cost and require optimization parameters tuning. As forour method, it is less complicated and more straightforward to control the size of the Voronoi cell. Though
the precision is slightly inferior, it signiﬁcantly improves
packing density.
3.3

Requirements of multi-sized particle sampling

Using the standard gridding approach (Fig.2) as an
example, when we generate particles in a unit cell, theminimum porosity is a ﬁxed value of 1 − π4 . It makes
minimum porosity become a constant value and impossible to adjust dynamically.
According to the properties of the power diagram, theplane are divided into several Voronoi cells and each ofthem will be a convex polygon with a weight factor that
can adjust its associated region. What inspires us is that
it is possible to obtain sampling results with diﬀerent
porosity sizes by adjusting each Voronoi cell’s weighting
coeﬃcients. During the Lloyd iteration, as changes in theweight coeﬃcients of each Voronoi cell impact the shapeand size of its associated region, it makes both the size ofIIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

particles and the pore size of its associated region changeable in each Voronoi cell. This property of the power
diagram makes it possible to achieve multi-sized particle
sampling within a given boundary. In order to achieve
multi-sized particle sampling with a higher packing density, i.e., ﬁnding a minimum porosity sampling case, wesummarize the requirements for the power diagram where
the minimum porosity is obtainable, which we consider
an ideal case:Requirement 1: If each side of the Voronoi cell has
the same length, then it is a regular polygon.
Requirement 2: If n Voronoi sites encircle a Voronoi
site, its associated Voronoi cell is a regular n sided polygon.
Requirement 3: In the case of satisfying n → ∞ andeach side d of Voronoi cell has the same length, the shapeof Voronoi cell is closer to the circle, porosity ϕ → 0 .
Requirement 4: Each Voronoi site must be at the centroid of its associated Voronoi cell.
Notably, the ideal case described above is only a scenario envisioned for each Voronoi cell. It is impossiblethat the Voronoi cell can fully satisfy all the given four
requirements locally or globally, but we can design theloss functions based on the requirements of an ideal case.
Speciﬁcally, as the requirements are formulated for each
Voronoi cell, the objective is to locally make each Voronoi
cell meet the requirements and ﬁnd a globally optimal solution.
3.4
3.4.1

Optimization
Objective function

To satisfy the requirements described in Section 3.3, we
sketched Fig.3 to represent an ideal case of the geometric
relationship between the Voronoi cells and their weights.
In this ideal case, Voronoi site A with B, C and D circling
it satiﬁes the conditions in Fig.4(b), which means thatcomputing the length of twin edge ab between Voronoi
sites A and B is suﬃcient to know the weight of anysite. Notably, it doesn’t mean that we ignore the weightof the other Voronoi site but it has been taken into account in the computing power diagram (the power diagram will be reconstructed before updating the weights
of Voronoi sites, as detailed in Algorithm 1). As for computing power diagrams, we mainly rely on the method

153

proposed by Aurenhammer34). They bring a general ddimentional power diagram up to compute with a d+1
dimensional convex hull. Brieﬂy, they transform a 2D
weighted Voronoi site s ∈ S into a sphere with coordi√
nates (xs , ys ) and radius r = ws mapped into the 3D
plane, and then intersect with other transformed Voronoi
sites on a 3D plane and ﬁnally apply an inverse transformation project onto the 2D plane (details can be foundin Aurenhammer34)or Nocaj et al.20)). The circles withthe solid line in Fig.3 are the projection results, and theirradii are the root of the weights of the associated Voronoi√
sites (i.e., the radius of Voronoi site A equal to wA ).
Assuming that there exist n points around Voronoi site
SA , Ideally, the side length dA of R (SA ) can be computed
as follows:qdA = 22
wA − r A(4)

where wA is the weight of Voronoi site SA and rA is theradius of Voronoi cell R (SA )’s maximum inscribed circle.
From the geometric properties of the power diagram,Equation 4 is only established in the case of Fig.4(b).
2, it will
When the weight of the Voronoi site is less than rAproduce the case (a) or case (c) in Fig.4, which means that
we will use a penalty function ψ (Section 3.4.2) insteadof shape loss function ΦP for preventing two sites frombeing too close or too far apart.
To satisfy requirements 1, 2, and 3 in Section 3.3, we
focus on information about the adjacent Voronoi sites
rather than the shape of the target Voronoi cell. It isbecause each edge of the target Voronoi cell is generatedusing the weights of the neighboring Voronoi sites. We
consider updating the weight of the target Voronoi site
for each iteration while the weights of the neighboring
Voronoi sites are ﬁxed. Consequently, to keep each edgeof the target Voronoi cell of the same length as possible,
the update tendency of the target Voronoi site weight canbe computed by incorporating the weights of all neighboring Voronoi sites. One can deﬁne the loss function that
controls the shape part of the porosity optimization as
follows:
ΦP (w, r) =X X

wi − wj + rj2 − ri2(5)i∈S j∈S−i

Meanwhile, we set a target particle radius distribution(Equation 10) in the initial stage to control the size of thesampled particles. According to this distribution function, we use the inverse transform method35)to ensurethat the generated target particle radii are uniformly dis-IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 4 Three possible cases between two neighboring Voronoi sites

tributed in each subinterval. Then we assign the generated target radius ri to each Voronoi site si when generating the initial power diagram or dynamically appending
particles (Section 3.4.3). Locally, the radius of the maximum inscribed circle in each Voronoi cell will approach
the user’s pre-set target radius. The global particle radius
distribution obtained from multi-sized particle sampling
will approximate the pre-deﬁned particle radius distribution. Therefore, the area of each Voronoi cell in the ideal
case (assume a regular n-sided polygon) can be computedby:Aideal
=
i

nri2
π
tan n

(6)

We re-generate the power diagram with new weights
with each iterative computation, and each Voronoi site
corresponds to a new Voronoi cell. In order to achieve
the target area for each newly generated Voronoi cell, wedeﬁne the loss function that controls the area part of theporosity optimization as follows:ΦA (A) =X Aideal− Acurrentiii∈S

Atotal(7)

where Acurrent
is the area of Voronoi site Si associiated region in the current iteration step, and Atotal isthe area of the boundary polygon Ω, which is the sumof all Voronoi cell areas. Noticeably, the Acurrent
is noti
modiﬁed directly by the cost function, but it depends on
the shape of the corresponding Voronoi cell in the updated power diagram. In addition, Acurrent
is also usedifor computing the porosity when we append new particles.
With joint the shape loss term ΦP and area loss termΦA , the ﬁnal loss function we minimise can be deﬁned as:Φ(w, r, A) = αΦP (w, r) + βΦA (A)power diagram to make each Voronoi cell meet the area
requirement. Moreover, as the geometric relationship between the Voronoi cells and their weights is known, it
allows us to optimize the shape of the Voronoi cells by
adjusting the weights. Balzer and Deussen ﬁrst proposed
optimally meeting area requirements via Lloyd’s method
19). It was followed by “Voronoi Treemaps”,20)which also
proved the practicability of this method. The core idea
is that the weight of each Voronoi site is increased or decreased in proportion to the missing or redundant areas
during Lloyd iteration. In our experiments, we found that
this method is also applicable for controlling the shapeof Voronoi cells, as shown in Fig.5 for the optimizationeﬀect of porosity. It should note that we could perform LBFGS method to optimize the power diagram32),33), but
Lloyd’s method is more straightforward, faster, and easier to implement. We would like to compare the impact
of Lloyd’s and L-BFGS methods on our experimental results in the future.
3.4.2 Constrained condition
While the Voronoi diagram ensures that each Voronoi
cell only contains one Voronoi site, the power diagram
may produce 3 cases, as shown in Fig.4, which showsthat the Voronoi cell may contain zero or several Voronoi
sites. Furthermore, the case of a Voronoi cell containing
no Voronoi site or more than one Voronoi site will leadto the wrong results of our multi-sized particle sampling
method. Our solution is not to perform the optimization
described in Section 3.4.1 when the distance dij between√√
Voronoi site i and j is greater than wi + wj or less√√than wi − wj . Instead, we deﬁne a penalty function
ψ for preventing two sites from being too close or too farapart, shown as follows:(8)

where α and β are two weighted coeﬃcients for adjusting global optimization. We apply Lloyd iteration to the154ψ (wi , dij , σk ) = σk d2ij − wi(9)where σk = k 2 presents the penalty factor in k-th iter-

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 5 The variation of porosity and loss value in the Bunny example

ation.
3.4.3

Dynamic particles appending

A fundamental issue with power diagram-based multisized particle sampling is that it requires searching forthe optimal number of particles to ﬁll the limited space.
The reason is that generating a power diagram is needed
to predetermine the number and weights of Voronoi sites.
While the optimum number of particles for ﬁlling the entire space can be roughly estimated from the pre-deﬁned
particle distribution, the weights of each Voronoi site
are diﬃcult to compute in advance. To solve this problem, we propose a dynamic particle appending based on
porosity. Speciﬁcally, we re-compute the porosity in thegiven boundary Ω at the end of each iteration and record
it. Within the recorded data series, if the current iteration step is Ti , we will extract the porosity data change.
Within a range from Ti−k to Ti (data length is k), we
determine whether to append particles. We are using theleast-squares method to ﬁt the extracted porosity dataand evaluate the slope of the ﬁtted line. When it is lower
than the threshold ϵ, we will append new particles according to the pre-deﬁned particle radius distribution. We set
k = 20 and ϵ = 1e−6 in all our experiments and illustrate the process of dynamic particles appending methodin Fig.6.

4.

Implementation

This section will introduce the multi-sized particle sampling method described in Algorithm 1. The input to ouralgorithm is a boundary polygon Ω, composed of a seriesof 2D vertices and a list of target particle radii, and thelist is generated from a user-deﬁned particle radius dis-155

tribution. The output is a power diagram consisting ofsampled n Voronoi sites, where each site contains a particle radius. First, before the porosity optimization, we use
the Lloyd iteration to ensure that each Voronoi site is located at the centroid of the Voronoi cell (line 22) followed
by a search for neighboring Voronoi sites (line 4). We do
not have to implement a speciﬁc neighbor search algorithm. The neighbor sites’ information can be obtained
from the QuickHull algorithm.
To address the 3 cases in the power diagram (Fig.4),
the weights of the target and neighboring Voronoi sites
determine whether to perform the objective function for
porosity optimization or a penalty function (lines 6 to
12). Line 13 describes the degree to which the area of
each Voronoi cell needs to be altered. Similarly, line 14
describes the degree of variation which is capable of making the shape of each Voronoi cell a regular n-sided polygon. Then we compute the values of the area and shapeweights respectively (weights wiA and wiP are increased or
decreased proportionally to ΦiA and ΦiP )for each Voronoi
cell that need to be updated for reaching the ideal case(lines 16 to 20). After updating the weighted factors, weneed to recompute the power diagram and perform Lloyd
iteration to make sure that each Voronoi site is locatedat the centroid of Voronoi cell (line 21 and line 22). Finally, we use “Dynamic Particles Appending Mechanism”
to determine the convergence of the current porosity andwhether the additional particles are needed (line 23). In
all of our experiments, we set α = 1 and β = 1 for astable and eﬀective optimization.
We have implemented our multi-sized particle samplingalgorithm using C++ on AMD Ryzen 9 3900X CPU. InIIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Algorithm 1 Multi-Sized Particle Sampling Based on
Porosity Optimization
Input: Boundary Polygon Ω, Target Particle Radius
R1 , R 2 , · · · , R n
Output: Power Diagram R(S), Voronoi Sites si ∈ S with aparticle radius
1: repeat
2:
R(S) ← Compute Power Diagram
3:
for each si in S do
4:
si Neighbor Searching(ni : Number of Neighbor Sites
sj ∈ Si′ ）5:ΦiP , Φtotal←0P
6:
for each sj √in Si′ do√√√7:
if (dij < | wi − wj |) ∨ ( wi + wj < dij ) theni28:ΦP += σk (dij − wi )9:else
P P (w −r2 )9:ΦiP += j ( j nji j − (wj − rj2 ))10:
end if11:
Φtotal+= ΦiPP
12:
end forAidea −Acurrent
n r2
i
13:
ΦiA = i
(Aidea= tani iπ )iniAtotal
14:

ΦiP =

ΦiP

ΦtotalP

15:
end for16:
for each si in S do
17:
wiA = 1 − (ΦiA − 1)218:
wiP = 1 − (ΦiP − 1)219:
wi += αwiA + βwiP
20:
end for
21:
R(S) ← Compute Power Diagram
22:
Performing Lloyd iteration
23:
Performing Dynamic Particles Appending Mechanism
24: until Φglobal < εFig. 6 Dynamic particles appending mechanism

line 2 of Algorithm 1, we use the quick hull algorithm36)
to compute the power diagram, which is stable, fast, andscalable for arbitrary dimensions. Additionally, we generate the straight skeleton37)in each Voronoi cell and compute the position and radius of the maximum inscribedcircle in the Voronoi cell based on the straight skeleton
as the output of Algorithm 1.

5.

Results

In the experiment, we set up the particle radius distribution at initialization by using a cumulative distributionfunction as follows:
0.5
Fradius (r) = 0.4
0.1

if 20 < r ≤ 30
if 30 < r ≤ 80

(10)

if 80 < r ≤ 150

At initialization, we randomly generate 100 points
within the boundary polygon Ω and set the maximum
iterations to 3000.
Figure 7 shows the results of our multi-sized particle sampling method at the initial stage, 200th iterations, 400th iterations, 2000th iterations, and 3000th iterations. We also recorded the values of the global porosity

156(Fig.5(a)) and the loss value (Fig.5(b)) during 3000 iterations by taking Bunny as an example. It is worth noting here that in Fig.5(b), each time the loss value drops
drastically, it is caused by triggering the dynamic particle
appending mechanism.
Even though the porosity in Fig.5(a) also declines
steadily when dynamic particle appending mechanism occurs, a sudden porosity increase occurs in the shortest
time after each particle appending. This phenomenon
is caused by the global porosity optimization having
reached a steady state before adding the particles. Nevertheless, after appending new particles, since the shape
optimization of each newly added Voronoi site is not sufﬁcient, this will result in the global porosity sudden increase. Additionally, we measured the predeﬁned particle
size distribution and optimized particle size distribution
over 3000 iterations in Fig.8. Intuitively, our optimized
particle size distribution is similar to the predeﬁned particle size distribution. Quantitatively, the error rate between optimized and predeﬁned particle size distribution
is roughly 15% according to RMSPE (Root Mean Squared
Percentage Error).
Except for the distributions we used in Fig.7, we also
experimented on whether it is possible to ﬁnd a minor
porosity under diﬀerent pre-deﬁned particle radius distributions. The top of Fig.9 presents the multi-sizedparticle sampling results after 3000 iterations, while thebottom part corresponds to the pre-deﬁned particle distribution. As a result of this experiment, setting a large

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 7 Results of our multi-sized particle sampling method

number of the small size particles will be helpful to obtaina more negligible porosity.

157

As the right-most case in Fig.9 shows that our method
allows to achieve a more complicated version of the parti-

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

cle distribution (containing three peaks), but in this caseit will be hard to promote the packing density.

6.

Discussion

Compared to uniform particle sampling method.
Both standard gridding approach5)and Poisson disk sampling4)can keep no overlapping between particles. Nevertheless, neither of them can improve the packing density
in a limited space while ensuring that the particles do not
overlap with each other. In addition, as these two methods do not obtain a balanced representation at the boundary, the sampled particle clouds hardly preserve the polygon’s original shape. Figure 7 demonstrates that the particles sampled by our method can ensure a smooth boundary and that there is no overlap between particles while
increasing the packing density (i.e., the global porosity
is eﬀectively reduced, as shown in Fig.5(a)). SPH-based
particle sampling methods6),7)can provide a uniform distribution at the boundary and produce a dense cloud ofparticles. Unfortunately, due to the specialty of the SPHbased method, it is not ensured that there is no overlapping between the sampled particles. It could be applied
to SPH-related algorithms, but it will produce unstable
simulation results when applied to DEM-related methods. Compared to the SPH-based method7), our algorithm could provide a stable DEM simulation, as shownin Fig.10.
Compared to multi-sized particle sampling
method. Most of the adaptive particle sampling methods1)–3)focused on the appearance of detailed ﬂuid surfaces or water splashes, and they had the property inwhich the particle radius became larger and closer to thecenter. As these methods are also based on SPH, they
have the same issues as SPH-based blue noise sampling7)

Fig. 8 Optimized and predetermined particle size
distribution

158. Furthermore, we also provide a comparison experimentwith the Poisson disk-based multi-sized particle samplingin Fig.10, which fails to achieve stable simulation results
due to the overlapping particle problem.
Compared to porosity with diﬀerent sampling
method. While we illustrate the instability simulation
when the Poisson disk sampling and SPH-based blue
noise sampling methods are applied to DEM in Fig.10,
the instability problem can be simply solved by detecting
the overlapping parts between particles and modifying
the particle radius. We admit that this approach is feasible but will induce a higher porosity, and the particle
radius distribution is hard to control. Nevertheless, ourapproach ensures no overlap between particles and promotes the packing density (minor porosity). To provethat our method can produce a relatively more minor
porosity, we adjust the size of each particle (the radiusof each particle is reduced by 1/2 the amount of overlap)
obtained by Poisson disk sampling and SPH-based blue
noise sampling so that there is no overlap between particles. As shown in Table 1, we statistically measure theporosity of 11 diﬀerent shapes of polygons respectively
by using Poisson disk sampling, SPH-based blue noise
sampling, and our method. Based on the results, the average porosity obtained by our sampling method is themost minor (0.165) and much better than the methods
Poisson disk sampling (0.697) and SPH-based blue noise
sampling (0.629).

7.

Conclusion

This paper proposed to improve the packing density ofparticle sampling in a restricted space as well as keeping no overlap between particles. Our method allows
users applying stable physics simulations based on DEM
methods and related frameworks. We have introduced amulti-sized particle sampling method based on porosity
optimization, which can obtain a particle distribution byminimizing the porosity in a pre-deﬁned 2D polygon. The
comparative experiments with alternative particle sampling methods demonstrate that we can deliver a stable
DEM-based sand simulation.
While the experimental results have proven that our
method can adapt to a variety of 2D polygons for multisized particle sampling, it still exists the case of undesirable sampling results as shown in Fig.11. These cases
showed that if a particle with a large radius preset neara narrow area, it will make the eﬀectiveness of local op-

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 9 Optimization results obtained from diﬀerent pre-deﬁned particle radius
distributions

timization quite unfavorable. Moreover, while our approach is designed to meet the pre-deﬁned particle radius
distribution by adjusting Voronoi cells shape and size continuously, it lacks how to keep the balance of the global
particle distribution, especially for the larger-sized particles. In the future, we would like to handle these problems
through an extra constraint on the particle position distribution. Furthermore, we expect to explore a way toextend our method to 3D multi-sized particle sampling
and adapt it to more physical simulation scenarios.

Acknowledgements
We would like to thank the anonymous reviewers
for their valuable suggestions and comments. This
work was supported by JSPS KAKENHI Grant Number JP20K11839 and JST SPRING Grant Number JPMJSP2124.

Fig. 10 Comparative experiments in sand simulation

References
1) B. Adams, M. Pauly, R. Keiser, L. J. Guibas:“Adaptively Sampled Particle Fluids”, ACM Trans. on Graphics, Vol.26, No.3,
pp.48:1-7 (2010).
2) R. Winchenbach, H. Hochstetter, A. Kolb:“Inﬁnite Continuous
Adaptivity for Incompressible SPH”, ACM Trans. on Graphics,
Vol.36, No.4, pp.102:1-10 (2017).
3) R. Winchenbach, A. Kolb:“Optimized Reﬁnement for Spatially Adaptive SPH”, ACM Trans. on Graphics, Vol.40, No.1,
pp.8:1-15 (2021).
4) R. Bridson:“Fast Poisson Disk Sampling in Arbitrary Dimensions”, Proc. of ACM SIGGRAPH Sketches, pp.22-es (2007).
5) M. Schwarz, H. P. Seidel:“Fast Parallel Surface and Solid Voxelization on GPUs”, ACM Trans. on Graphics, Vol.29, No.6,
pp.179:1-10 (2010).

1596) H. Schechter, R. Bridson:“Ghost SPH for Animating Water”,
ACM Trans. on Graphics, Vol.31, No.4, pp.61:1-8 (2012).
7) M. Jiang, Y. Zhou, R. Wang, R. Southern, J. -J. Zhang:“Blue
Noise Sampling Using an SPH-based Method”, ACM Trans.
on Graphics, Vol.34, No.6, pp.211:1-11 (2015).
8) F. De Goes, C. Wallez, J. Huang, D. Pavlov, M. Desbrun:“Power Particles: An Incompressible Fluid Solver Basedon Power Diagrams”, ACM Trans. on Graphics, Vol.34, No.4,
pp.50:1-11 (2015).
9) X. Zhai, F. Hou, H. Qin, A. Hao:“Fluid Simulation with Adaptive Staggered Power Particles on GPUs”, IEEE Trans. on
Visualization and Computer Graphics, Vol.26, No.6, pp.22342246 (2018).
10) M. Aanjaneya, M. Gao, H. Liu, C. Batty, E. Sifakis:“Power Di-

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Table 1 Table of porosity data with diﬀerent particle sampling method
Example Name Poisson Disk Sampling SPH Based Blue Noise Sampling Ours
Dragon
0.677
0.616
0.18
Horse
0.625
0.602
0.16
Lucy
0.605
0.613
0.172
Woody
0.703
0.598
0.161
Homer
0.708
0.655
0.172
Cheburashka
0.706
0.599
0.17
Teapot
0.713
0.635
0.164
Spot
0.772
0.653
0.168
Cow
0.704
0.649
0.17
Armadillo
0.695
0.604
0.174
Bunny
0.761
0.704
0.124
Average
0.697
0.629
0.165

Fig. 11 Some of the failure cases
agrams and Sparse Paged Grids for High Resolution Adaptive
Liquids”, ACM Trans. on Graphics, Vol.36, No.4, pp.140:1-12(2017).

20) A. Nocaj, U. Brandes:“Computing Voronoi Treemaps: Faster,
Simpler, and Resolution-independent”, Computer Graphics
Forum, Vol.31, No.3, pp.855-864 (2012).

11) N. Ray, D. Sokolov, S. Lefebvre, B. Lévy:“Meshless Voronoi on
the GPU”, ACM Trans. on Graphics, Vol.37, No.6, pp.265:1-12(2018).

21) W. Rungjiratananon, Z. Szego, Y. Kanamori, T. Nishita:“Realtime Animation of Sand-water Interaction”, Computer Graphics Forum, Vol.27, No.7, pp.1887-1893 (2008).

12) X. Liu, L. Ma, J. Guo, D. -M. Yan:“Parallel Computation of
3D Clipped Voronoi Diagrams”, IEEE Trans. on Visualization
and Computer Graphics, Vol.28, No.2, pp.1363-1372 (2020).

22) T. Lenaerts, B. Adams, P. Dutré:“Porous Flow in Particlebased Fluid Simulations”, ACM Trans. on Graphics, Vol.27,
No.3, pp.49:1-8 (2008).

13) J. Zheng, T. -S. Tan:“Computing Centroidal Voronoi Tessellation Using the GPU”, Proc. of the ACM Symposium on Interactive 3D Graphics and Games, pp.14:1-9 (2020).

23) T. Lenaerts, P.Dutré:“Mixing Fluids and Granular Materials”,
Computer Graphics Forum, Vol.28, No.2, pp.213-218 (2009).

14) J. Basselin, L. Alonso, N. Ray, D. Sokolov, S. Lefebvre, B.
Lévy:“Restricted Power Diagrams on the GPU”, Computer
Graphics Forum, Vol.40, No.2, pp.1-12 (2021).
15) M. Balzer, T. Schlömer, O. Deussen:“Capacity-Constrained
Point Distributions: A Variant of Lloyd’s Method”, ACM
Trans. on Graphics, Vol.28, No.3, pp.86:1-8 (2009).
16) M. Balzer:“Capacity-Constrained Voronoi Diagrams in Continuous Spaces”, Proc. of the IEEE Symposium on Voronoi
Diagrams, pp.79-88 (2009).
17) Y. Liu, W. Wang, B. Lévy, F. Sun, D. -M. Yan, L. Lu, C.
Yang:“On Centroidal Voronoi Tessellation―Energy Smoothness and Fast Computation”, ACM Trans. on Graphics, Vol.28,
No.4, pp.101:1-17 (2009).
18) R. Chen, C. Gotsman:“Parallel Blue-noise Sampling by Constrained Farthest Point Optimization”, Computer Graphics Forum, Vol.31, No.5, pp.1775-1785 (2012).
19) M. Balzer, O. Deussen, C. Lewerentz:“Voronoi Treemaps for
the Visualization of Software Metrics”, Proc. of the ACM Symposium on Software Visualization, pp.165-172 (2005).

16024) R. Narain, A. Golas, M. -C. Lin:“Free-ﬂowing Granular Materials with Two-way Solid Coupling , ACM Trans. on Graphics,
Vol.29, No.6, pp.173:1-10 (2010).
25) G. Klár, T. Gast, A. Pradhana, C. Fu, C. Schroeder, C.
Jiang, J. Teran:“Drucker-prager Elastoplasticity for Sand Animation”, ACM Trans. on Graphics, Vol.35, No.4, pp.103:1-12(2016).
26) T. Yang, J. Chang, B. Ren, M. -C. Lin, J. -J. Zhang, S. M. Hu:“Fast Multiple-ﬂuid Simulation Using Helmholtz Free
Energy”, ACM Trans. on Graphics, Vol.34, No.6, pp.201:1-11(2015).
27) T. Yang, J. Chang, M. -C. Lin, R. R. Martin, J. -J. Zhang, S. M. Hu:“A Uniﬁed Particle System Framework for Multi-phase,
Multi-material Visual Simulations”, ACM Trans. on Graphics,
Vol.36, No.6, pp.224:1-13 (2017).
28) A. P. Tampubolon, T. Gast, G. Klár, C. Fu, J. Teran, C.
Jiang, K. Museth:“Multi-species Simulation of Porous Sand
and Water Mixtures”, ACM Trans. on Graphics, Vol.36, No.4,
pp.105:1-11 (2017).
29) M. Gao, A. Pradhana, X. Han, Q. Guo, G. Kot, E. Sifakis,

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）
C. Jiang:“Animating Fluid Sediment Mixture in Particle-laden
Flows”, ACM Trans. on Graphics, Vol.37, No.4, pp.149:1-11(2018).
30) Y. Yue, B. Smith, P. -Y. Chen, M. Chantharayukhonthorn,
K. Kamrin, E. Grinspun:“Hybrid Grains: Adaptive Couplingof Discrete and Continuum Simulations of Granular Media”,
ACM Trans. on Graphics, Vol.37, No.6, pp.283:1-19 (2018).
31) X. Wang, M. Fujisawa, M. Mikawa:“Visual Simulation of SoilStructure Destruction with Seepage Flows”, Proc. of the ACM
on Computer Graphics and Interactive Techniques, Vol.4, No.3,
pp.41:1-18 (2021).
32) L. Zheng, Y. Yao, W. Wu, B. Xu, G. Zhang:“A Novel Computation Method of Hybrid Capacity Constrained Centroidal
Power Diagram”, Computers and Graphics, Vol.97, pp.108-116(2021).
33) S. -Q. Xin, B. Lévy, Z. Chen, L. Chu, Y. Yu, C. Tu, W.
Wang:“Centroidal Power Diagrams with Capacity Constraints:
Computation, Applications, and Extension”, ACM Trans. on
Graphics, Vol.35, No.6, pp.244:1-12 (2016).
34) F. Aurenhammer:“Power Diagrams: Properties, Algorithms
and Applications”, SIAM Journal on Computing, Vol.16, No.1,
pp.78-96 (1987).
35) J. M. Steele:Non-Uniform Random Variate Generation,
Springer (1987)
36) C. B. Barber, D. P. Dobkin, H. Huhdanpaa:“The Quickhull
Algorithm for Convex Hulls”, ACM Trans. on Mathematical
Software, Vol.22, No.4, pp.469-483 (1996).
37) P. Felkel, S. Obdrzalek:“Straight Skeleton Implementation”,
Proc. of Spring Conference on Computer Graphics, pp.210-218(1998).

(Received March 10, 2022)(Revised July 7, 2022)

Xu

WANG

He received his BS degree in information management and information system from Shanghai University in 2016, and MS degree in media science from Tokyo University of Technology in 2020. He is currently a PhD candidate at Physics Based Computer Graphics
Laboratory, University of Tsukuba. His research interests include physics-based simulation, style transfer, and all the relevant topics
in computer graphics.

Makoto FUJISAWA (Member)
He is currently an Associate Professor in
the Faculty of Library, Information and Media Science, University of Tsukuba since2021.
He received B.Eng., M.Eng., andPh.D. degrees in mechanical engineering from
Shizuoka University in 2003, 2005, and 2008
respectively. He worked for Nara Instituteof Science and Technology from 2008 to 2010
and University of Tsukuba from 2011 to 2020
as an assistant professor. His research interests include computer graphics and physics
simulation. He is a member of ACM, IEEE
CS, IIEEJ, IPSJ and VRSJ.

161

Masahiko MIKAWA
He is currently an Associate Professor in the
Faculty of Library, Information and Media
Science, University of Tsukuba, Japan since2006. He received B.Eng., M.Eng., and Ph.D.
degrees from Osaka University in 1992, 1994
and 2001 respectively. He worked for NTT
Access Network Systems Laboratories from1994 to 2001, NTT Service Integration Laboratories from 2001 to 2003 and was a Lecturer in the Graduate School of Library, Information and Media Studies, University ofTsukuba from 2003 to 2006. He is a member
of RSJ, SICE, SOFT and IEEE.

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）-- Special Issue on Visual Computing --Real-Time Rendering of Oil Film with Flexible Properties
Subroto Prasetyo HUDIONO†, Tomoya ITO†† (Member), Yuriko TAKESHIMA† (Member), Tsukasa KIKUCHI† (Member)†Tokyo University of Technology, ††Hachinohe Institute of Technology<Summary> Thin-film interference is widely known for the reason of the iridescence on a bubble and spilled oil. Although
it is not difficult to calculate the interference of light waves, calculating high precision interference is still considered
impractical in real-time rendering. The proposed method is to separate the interference calculation by pre-computing
calculation and save the value in lookup texture. But instead of using angle and thickness like other methods, we are using
optical path difference as an index in lookup texture to create a flexible and compact lookup texture. The rest of the calculation,
which calculates the optical path difference from the thin-film attribute will be done in real-time. Compared with other methods,
the lookup texture generated using optical path difference will consume significantly less memory since it reduces the texture
dimension to a one-dimensional texture. With this, we can handle different parameter changes such as refractive index andreflectance value without generating a new lookup texture or impacting the rendering performance.
Keywords: real-time rendering, thin-film interference, iridescence

is often seen in thin-film interference. Our approach relieson this pattern to avoid recalculating the color in every
frame. By considering the physical parameters, we focus on
the pre-computation of the cause of interference to make it
flexible when the parameters such as refractive index, film
thickness, and reflectance value change. With this, we candesign an efficient look-up texture that consumes lessmemory with a wide range of usage.
In this paper, we will discuss the calculation that isneeded to render a high accuracy interference and theimportance of each thin-film parameter. From there we canseparate the calculation between precomputed interference
calculation and real-time attribute calculation. In addition, to
produce realistic results, the rendering process from thecreation of lookup texture to the real-time rending is also
explained in detail. Finally, we will compare the results inmemory consumption and calculation speed to get the
improvement of our method from the previous methods.

1. Introduction
Iridescent surfaces are quite common in nature and has abeautiful property where the surfaces change their colors
depending on the direction from which they are viewed.
This occurs because light wave interference happens at thesurface and causes some spectra to dim or even disappear.
The most common cause is thin-film interference, a lightwave interference caused by the thin film around the objects.
We often see such phenomenon on soap bubbles or spilledoil, where a gradient of color appears on the surface if seenfrom different angles.
In computer graphics, iridescent surfaces are also often
used as a material of an object, either to enhance the beautyof the objects or make them seem more realistic. It is often
used as a material coating for cars or glasses. These
examples can be rendered quickly since usually they have afixed structure. Other common examples are oil spills andsoap bubbles. Because of their liquid properties, thethickness of the thin film keeps changing due to surface
tension. Different types of oil also have different refractive
indexes, affecting the resulting color of the surface. These
properties make the rendering too complex to be
pre-computed while rendering it in real-time is too
expensive and inefficient.
While iridescence seems to have a wide range of colors,
the colors produced are due to the missing spectra after theinterference. Therefore, a certain sequential pattern of color

2. Related WorkThin-film interference is often seen in glass coatings orsoap bubbles. Most studies use soap bubbles to observe thisbehavior, and many have tried to reproduce it in computer
graphics, with some excellent results attained by simulating
bubble structures1) and surface tension flow on the bubble
surface2). However, such studies focused more on the bubble
itself instead of the interference.
In the field of thin-film interference, one of the earliest

162

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

models in computer graphics was that of Smits and Meyer3).
Their model can accurately reproduce iridescent effects fromthin films of various thicknesses. However, their method islimited to certain conditions, for example, it only applies to asmooth surface, and it does not consider inter-reflections
inside the thin film. These limitations are handled splendidlyin the model of Belcour and Barla4), which uses microfacet
theory for rendering the iridescent effects of varying
thickness. Some works are related to iridescent color but use
different physical structures. For instance, the model of Icart
and Arquès5) divides the intensity of reflected light into
several intensities. A unique approach was proposed bySun6), whose model uses multilayer interference to simulate
natural patterns found in animals and insects such as on thewings of morpho butterflies.
Despite much research on thin-film interference, there
have been a few studies on its real-time rendering. One ofthem was an approach by Iwasaki et al.7), which could
address this problem by calculating the interference inadvance and saving it as a two-dimensional texture. This
pre-computation approach works well for soap bubbles
where the thickness of the thin film keeps changing.
However, it was inflexible since a new texture need to berendered when there are changes in the parameter such as
refractive index. Moreover, there are many duplication
calculations due to using film thickness and incidence angle
as a texture index. This makes the texture for interference
needs too much size to render high accuracy colors.
To solve this problem, we developed a new model thatcan render iridescence more accurately while using lessmemory. Our approach focuses on creating an interference
map based on the cause of the interference, which is thedifference in the optical path, instead of the incidence angle
and film thickness. With this approach, we can remove theinfluence of film properties and the angle of incoming lightfrom the interference map. This makes it more flexible whenrendering multiple objects with different film properties.

Fig.1

Light paths inside an oil film

different distances, the wavelength may be shifted, causing
some frequencies to interfere and become dimmed or
disappear, and leading to the reflected light missing somecolor spectra. This phenomenon is called destructive
interference; constructive interference is when both waves
having the same optical path merge and cause a brighter
color.
3.1 Interference between two waves
Interference between light waves happens because thewave inside the thin film travels a different distance andcauses shifting when merged with other waves on thesurface of the film. This difference in distance is called
optical path difference (OPD). Film thickness and lightangle are the main properties that determine OPD. Depends
on the direction of the light comes, the light travels farther
inside the film. This is the reason why we see different
colors when we see a thin-film surface from a different
angle.
The refractive index is also an important parameter inthin-film interference. It is a parameter to determine thespeed of light when traveling inside a substance. For
example, light travels at 1.5 slower inside an oil film with arefractive index of 1.5. The difference in refractive index
between substances also causes the light to change direction
or refract. The relationship between incidence angle,
refractive angle, and refractive index is defined in Snell’slaw,

3. Physical Model of Thin Film Interference

sin θ
sin θ

The low surface tension of oil and its hydrophobic
properties make it spread easily on objects or water surfaces
under the influence of gravity. It can spread into a thin filmthat causes light waves to interfere with each other. Thisphenomenon is called thin-film interference, whereby lightwaves are reflected by the upper layer interfere with thediffracted light reflected by the bottom layer of the thin film,
as illustrated in Fig. 1.
Since the light reflected by the bottom surface travels

n
,
n

1

where θ and n are the incidence angle or angle from thelight source and the refractive index of the first medium, andθ and n are the refractive angle and the refractive index
of the film.
The Optical Path Difference or OPD is written as λ,because we need to find the phase difference, the distancedifference between wave needs to be measured with the163IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

wavelength. With λ as the wavelength, andrepresents
the number of waves difference. The OPD can be foundwith the formula
2

cos θ

λ,

2

where
is the thickness of the thin film. From the OPD we
can determine whether the interference is constructive ordestructive. If
is an integer, the difference is a multipleof λ, and therefore the waves inside and outside of the film
are aligned and the interference is constructive. Meanwhile,if
has a value of 0.5, the wave shifts by half of its length,
and the interference is destructive, resulting in the wave
disappearing.
Phase shift is also a common occurrence in thin-film
interference. When light is reflected at the surface of the film,
if the refractive index of the film is lower than that of themedium through which the light passes, a phase shift will
occur. This results in a wave shift of 180° or half of the
wavelength. This also affects OPD, which increases by 0.5λ
if there is a phase shift.
With the OPD, we can calculate the intensity after
interference at the surface of the thin film. We can define anequation between reflected wave and refracted wave as,
x, λ

cos

cos

cos 0

cos

cos

λ !,

λ ,

λ ,

3

where
is the peak amplitude of the light wave before
entering inside the film,
is the reflectance,is thetransmittance value,is the wavenumber,is the initial
phase of the wave, and λ is the OPD. Since we use thedifference in the optical path, interference can be assumed at
the start of angular frequency which makes
to be 0.
In the actual case, the light that enters inside the film isreflected and attenuated several times until losing its energy.
Every time a portion of the light travels out of the film, it has
traveled a certain distance that depends on how many times
it was reflected inside the film. Thus, interference happens
multiple times at different phases as illustrated in Fig. 1.
With the assumption that the thin film is on a solid surface
that fully reflects the light, we can define it as the equation,
λ

cos

cos 3

λ

λ

2π
cos $
λ&
λ

cos $3
(

,∞

)

*-

2π
λ&
λ

*+

cos 2

⋯,

⋯

cos 2πwhere the Σ represents the total radiance from all interfered
waves andrepresents the number of lights that went
outside the film. Since the material color is defined byalbedo texture, the reflectance value of the material assumed
to be 1 and can be removed from the equation.
3.2 Visible spectrum
The visible light spectrum covers wavelengths ranging
from about 380 to 750 nanometers. Each wavelength
represents a different color, while white light is acombination of multiple lights of different wavelengths.
When the light gets inside the thin film, it contains a lot ofdifferent waves at different wavelengths. Therefore, we cansay that the interference happens in all the visible light
spectra that got inside the thin film.
We can accurately calculate the output color bydetermining the radiance of each light that enters inside thethin film. This can be defined as,
0

λ

λ

5

where λ is the spectral radiance of the light at a certain
wavelength. The number of samples that we use in thisinterference equation also affects the accuracy of the final
color. As shown in Fig. 2, fewer samples tend to result inover-saturated color output, thus a higher number of wave
samples is required to accurately render the color whencalculating the total radiance.

λ2π
cos $2λ&λFig.2 A difference in saturation between rendering the interference with 10
wave samples (a) and 370 wave samples (b)4. Rendering Method
Simulating a large number of wavelength samples to
produce realistic output colors in thin-film interference iscomputationally expensive. Moreover, rendering sRGB
color space using wavelengths with high accuracy could be
impractical in real-time rendering. Especially when there arevariations in viewing angle and thin-film properties, the

4

.

164

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig.5

Fig.3

An OPD texture (a) takes 640 bytes and a 2D texture (b) takes
4,096 bytes of memory

Fig.4

Two-dimensional interference map from thicknesses and incidence
angles without considering the reflectance value

output color will have to be recomputed4).
To date, the simplest way to render thin-film
interference in real-time rendering has been to pre-compute
all possible angles and thicknesses from thin-film
interference use them as a texture7). Although this could be
used in real-time rendering, it is not flexible enough tohandle the change in parameters such as refractive index andreflectance value. Moreover, it contains duplicate values
along the Y-axis and uses memory consumption
ineffectively. As seen in Fig. 3, 2D texture can consume
more memory with lower accuracy compared to using OPD
texture.

Lookup texture by calculating the interference with every possible
OPD and reflectance value

the texture will save multiple duplicate values. This can beseen in Fig. 4. By ignoring the reflectance value to make thepattern visible, the value is mostly the same between the Yaxis, only stretched because of the incidence angle.
In generating an interference map, we calculateinterference with the formula from Eq. (4) to simulate theinterference that happens between the light inside the film.
This interference between light waves has less effect withlonger OPD, as we can see in Fig. 4. From this, we cansimulate the interference that happens with OPD
below 5μm.
The input for the calculation will be the OPD λ andthe wavelength
λ . To simulate the interference
accurately, we sample 370 wavelengths from 350 to 750.
With the same OPD on every wavelength, we will get
different multiplication numbersthat determine
whether there is destructive interference or not. Finally,
we average the results after changing the wavelength to
RGB colors and save it in the texture.
4.2 Reflectance value
Another factor that can affect interference is thereflectance value of the thin film. Reflectance value is the
proportion of the wave when separated into a refracted wave
and a reflected wave when hitting the thin-film surface.
Since the reflectance value is the proportion of the divided
light, we can represent it as a scale number from 0 to 1.
With reflectance value, the lookup texture will become a2-dimensional texture with reflectance value as the Y index.
The result of this texture can be seen in Fig.5. We will havea bigger lookup texture, but since the reflectance value islinear and doesn’t need high precision, we can use a smallvalue with the incremental of 0.1 precision. Moreover, theresult is horizontally symmetrical, so we only need to
calculate half of the value which is 0 to 0.5. This will
increase the lookup texture by 6 times larger, but still muchsmaller by using a fully 2-dimensional texture.

4.1 Optical path difference
We propose a new rendering method in which, insteadof using thickness and angle to determine the resulting color
in a two-dimensional texture, we map the output color in a
one-dimensional texture and use OPD as an index.
The main cause of differences in color output inthin-film interference is the shift in wavelength due to thedifferent distances traveled, or OPD. First, the incidence
angle and film thickness are input to determine the OPD,
however, using both of them as an index is inefficient since165

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

4.3 Wavelength to RGB color space
The output of interference calculation (Eq. (4)) is theradiance of each wavelength, from this we use CIE's color
matching functions to transform the wavelength (λ) into
RGB8). This color-matching function is based on theequation,
0

λ

L

0

λ M λ λ,N0

λ O λ λ,K
K

Results of thin-film interference photographed for a real soap film
(left10) and middle1)) and our simulated oil film (right)

Fig.7

Thin-film interference implemented in a cube with different
materials

λ λ,

J

K

Fig.6

7

Given L λ is the spectral radiance of the wavelength λ.
This tristimulus value (XYZ) can be transformed to RGB
using the following matrix9),
Q
(R .
S

3.24048 T1.53715 T0.49854 J(T0.96926 1.87599
0.04156 . (L .,0.05565 T0.20404 1.05731 N

8Then the RGB value will be scaled by the number ofwavelength samples used to simulate the interference.

5. Results
The accuracy of our approach was tested by comparing
the result of interference color with a photograph of a realsoap film from Atkin and Eliot10) (Fig. 6). With color
temperature adjustment, our simulated oil film, rendered in a
linear incremental of film thickness, could well reproduce
the color pattern of the photograph of interference. Since ourapproach is using simulated optical path difference tocalculate the interference, the results would be the sameeven for materials with different refractive indexes.
We have implemented our approach in GLSL shaders
using our custom renderer running in Vulkan with AMD
Radeon Pro 450. We tested the speed performance of
calculating OPD in real-time rendering with physically
based rendering (PBR) on a metallic surface object. The
rendering times for each frame are given in Table 1.
The first row shows the rendering times without
interference for a normal object and material by PBR. The
second row shows the rendering time when the interference

was calculated directly while rendering the scene. The last
two rows show the results of testing the scene with 2D
format texture and our approach. From this table, we knowthat a naive approach to calculating the thin-film
interference with high accuracy isn’t feasible in real-timerendering. And while our approach takes more rendering
time the difference is too small and didn’t scale by thecomplexity of the scene
The memory consumption of OPD texture is alsosignificantly smaller compared with 2D texture lookup.
Because of this, we can increase the resolution withoutworrying about the size of the texture lookup. We calculate
the estimated size of the texture in Table 2

Table 1 Rendering time per frame in millisecond

Table 2 Memory consumption from different sample size

SphereBunny

TextureBubbleResolution

2D Texture

No Interference

8.8 ms12.8 ms

28.0 ms

15.2 ms

1024 samples

4 MB

4 KB

24 KB

OPD Texture OPD with ReflectiveNaïve

134.6 ms

211.8 ms

114.9 ms

147.8 ms

2048 samples

8 MB

8 KB

48 KB

Texture (2D)

9.6 ms

13.4 ms

28.3 ms

16.8 ms

8192 samples

32 MB

32 KB

192 KB

Ours (1D)

10.1 ms

13.6 ms

28.3 ms

17.0 ms

16384 samples

64 MB

64 KB

384 KB

166

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig.8

Fig.9

iron oxide will also increase, forming a thin film and causing
the steel color to change. Therefore, we tried to reproduce
the thin-film interference on the tempered steel by
increasing the thickness of the film from 0.12μm to
0.24μm.
To test the importance of refractive index as a parameterof thin films, we tried to compare the difference whenincreasing the refractive index and the film thickness in a
sphere. From the test in Fig.10, we know that the change infilm thickness has a different effect from the change inrefractive index.
Increasing the thickness of the film will give more
gradient and make the material more sensitive to color
change, this is because by increasing the thickness, the OPD
will also increase and ended up using a wider area in the
lookup texture. Meanwhile, the change in the refractive
index gives us shorter gradient colors. This is because therefractive index also affects the refractive angle and causes
the view angle to have a smaller effect on the OPD.
One of the examples of an object with a high refractive
index is tempered steel. Tempered still tends to have thesame color even though viewed from different angles. By
reducing the refractive index, the area of the lookup tablethat is used to get the color is also scaled down. This result
cannot be achieved when using a 2D lookup texture basedon film thickness and incidence angle.

Rendered soap bubble with thin-film interference

Tempered steel color reproduced by increasing the film thickness
to simulate the oxidation

6. Conclusion and Future Works

In 2D texture from Iwasaki et al. 7), the texture is
generated from the film’s thickness and incidence angle.
Since both are important factors for color accuracy, theyneed to have high accuracy for generating a lookup texture.
With this, we generate the lookup texture with at least 1024
samples for the incidence angle with accuracy per 0.1
degree.
We tested our method with physically based rendering
on different materials to test the color accuracy (Fig. 7).
While the surface that has a thin film will reflect light, someparts of the material have low reflectance, and might notshow any color differences. With that in mind, we tested ourapproach by implementing the interference on the part of theobject which reflects light.
We also reproduced objects that have a natural thin-film
interference property like soap bubble (Fig. 8) and tempered
steel (Fig. 9). The heating process in tempered steel cancause a layer of iron oxide to form on the surface of the steel.
As the temperature of the steel increases, the thickness of the167

We have presented a new approach to generate lookup
texture with optical path difference as an index. One of theadvantages of our method is the flexible properties of a thinfilm that can be changed without generating a new lookup
texture. This includes the film thickness, refractive index,
and reflective value, while the texture used from othermethods is only able to change the film thickness.
Moreover, with the pre-computation of the interference,
there will be no limitation in the first rendering time, and amore accurate RGB model could be used to achieve moreaccurate color when transforming wavelength into color
without affecting the rendering time. This could be usefulwhen rendering spilled oil on a road or car material in games
where the scene needs to have a high frame rate.
Our contribution mainly focuses on providing a new
flexible way to apply iridescent materials to different typesof objects with high accuracy in real-time rendering. As aresult, a scene with different types of objects with different
material properties could be rendered accurately by only
using one small interference texture in real-time. It is also

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig.10The difference in changing the thickness and refractive index of the thin film. The left image is the sphere with a high refractive index, themiddle image is the default parameters, and the right image is the image with a thicker thin film

Referencespossible to implement the proposed method directly into
production without any adjustments.
The limitation of our approach is the assumption that theobject only has a single film and a smooth surface in which
the surface will fully reflect the light. When the object has
multiple films, the interference becomes much more
complex that the interference map might not be possible toproduce an accurate result. One of the reliable methods is torender all the interference in real-time but it is not possible
for real-time rendering. But it is rare to find an object withmore than one thin film.
Future work includes rendering multiple layers of thinfilms. The proposed model focused mainly on rendering one
film on the object's surface. Meanwhile, in a spilled oil, thelight that passed through also interact with the water below
the oil film making the object has two films with different
refractive index and different thickness. This could be aninteresting topic to study considering spilled oil is a common
occurrence in the real world.
Finally, even though real-time rendering techniques have
improved greatly together with the development ofcomputer hardware in recent years, providing a real-timerendering with high realism is still a challenge. A lot of
aspects need to be considered, and even small improvements
could help to save up time during rendering. In the case ofthin-film interference, the color outputs already have great
accuracy, so instead of the output color, a film behavior likehow the thickness move in oil film due to the surface tension
is also an interesting part of future work.

1)

S. Ishida, P. Synak, F. Narita, T. Hachisuka, C. Wojtan: “A Model for
Soap Film Dynamics with Evolving Thickness,” ACM Trans. on
Graphics, Vol.39, No.4, pp. 31:1–31:11 (2020).

2)

W. Huang, J. Iseringhausen, T. Kneiphof, Z. Qu, C. Jiang, M. B. Hullin:“Chemomechanical Simulation of Soap Film Flow on Spherical
Bubbles,” ACM Trans. on Graphics, Vol.39, No.4 (2020).

3)

B. Smits, G. Meyer: “Newton’s Colors: Simulating Interference4)

L. Belcour, P. Barla: “A Practical Extension to Microfacet Theory for

Phenomena in Realistic Image Synthesis,” (1992).

the Modeling of Varying Iridescence,” ACM Trans. on Graphics,
Vol.36, No.4, p. 14 (2017).
5)

I. Icart, D. Arquès: “An Illumination Model for a System of Isotropic
Substrate- Isotropic Thin Film with Identical Rough Boundaries,” In
Proceedings of the 10th Eurographics conference on Rendering,
Eurographics Association, Goslar, DEU, 261–272 (1999).

6)

Yinlong Sun: “Rendering Biological Iridescences with RGB-Based
Renderers,” ACM Trans. on Graphics, Vol.25, No.1, pp. 100-129(2006).

7)

K. Iwasaki, K. Matsuzawa, T. Nishita: “Real-time rendering of soap
bubbles taking into account light interference,” Proceedings Computer
Graphics International, pp. 344-348 (2004).

8)

C. Wyman, Peter-Pike J. Sloan, P. Shirley: “Simple Analytic
Approximations to the CIE XYZ Color Matching Functions,” Journalof Computer Graphics Techniques, Vol.2, No.2 (2013).

9)

R. Durikovic, Ryou Kimura: “GPU Rendering of the Thin Film on
Paints with Full Spectrum,” Tenth International Conference on
Information Visualisation (IV’06), pp. 751-756 (2006).

10) L. J. Atkins, R. C. Elliott: “Investigating thin film interference with adigital camera,” American Journal of Physics, Vol.78 (2010).

(Received March 11, 2022)

168

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Subroto Prasetyo HUDIONOHe received his B. Sc. degree from Institut Sains
dan Teknologi Terpadu Surabaya in 2014. After
entered Tokyo University of Technology,
Department of Media Science in 2020, hecompleted his master’s degree program in 2022. His
research interests include computer graphics, game,and real-time rendering.

Tomoya ITO(Member)

He received the Ph.D degree from Iwate
University, Japan in Engineering. He has been a
Professor in the Department of Engineering, atHachinohe Institute of Technology since 2004. His
research interests include computer graphics, game,and media art. He is a member of Information
Processing Society of Japan, Game Amusement
Society, and the Society of Art and Science.

Yuriko TAKESHIMA(Member)

She received her B. Sc., M. Sc., and Ph.D.
degrees in Computer Sciences from Ochanomizu
University in 1994, 1996, and 1999, respectively.
She is currently a Professor at School of Media
Science, Tokyo University of Technology. Her
research interests include volume visualization,
topology-based visualization, and visualization
environments.

Tsukasa KIKUCHI(Member)

He received his Ph.D degree from Iwate
University in Engineering in 1999. He is currently aProfessor at School of Media Science, Tokyo
University of Technology. His research interests
include computer graphics, visual simulation,
procedural animation and contents design. He is amember of Information Processing Society ofJapan, and the Society of Art and Science.

169

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）
IIEEJ

Paper
Eyeliner Looked “Attractive” by Eyelids Shape Using Eyeliner Models
Rena OKURI† , Shuhei KODAMA† , Tokiichiro TAKAHASHI† †† (Fellow )† Tokyo Denki University , †† ASTRODESIGN,Inc.

<Summary> “Eyeliner,” a makeup method makes the eyes look larger, and makes the face so-called‘attractive.’ Eyelids shapes of Asian women are roughly categorized into “double,” “hooded,” and “single.”
We assumed that the most attractive eyeliner for her would depend on her eyelids shape. In this paper, wequantitatively measured the length and the thickness of eyeliner (eyeliner style) that make the eyes appear
the largest without any sense of incongruity for each eyelids shape and investigated whether the eyeliner
style is attractive. The purpose of this paper is to clarify attractive eyeliner styles by eyelids shape. The
measurement was conducted using digital makeup with eyeliner models devised in this study. The results
showed that the eyeliner that makes the eyes look the largest without any sense of incongruity is the most
attractive eyeliner for each eyelids shape.
Keywords: makeup, eyeliner, eyeliner model

1.

Introduction

Facial impression is one of the important factors in impressing personal so-called ‘attractiveness.’ According toDipboye et al.1), an attractive person even has advantages
in social life. Facial features are related to facial attractiveness and faces with features of high symmetry andaverageness are judged to be attractive2). Furthermore,according to experimental studies on facial morphology,
large eyes have one of the features that increases the attractiveness of a woman face2) because they give the impression to be feminine and youthful3),4).
It is diﬃcult to change the shape of the face itself tohave these features. However, makeup, applying various
cosmetics to the face, is a way to approach the face withattractive features only in appearance without changing
the shape of the face. Makeup has the eﬀect of increasing
facial attractiveness and beauty5)–7).
There have been developed various makeup methods,
and examining their eﬀectiveness is an important role toimprove the attractiveness with the recent growing interest in facial beauty and attractiveness. Previous work on
the eﬀects of makeup have shown that eye makeup is particularly eﬀective in increasing attractiveness by changing
the apparent features of the eyes8).
Eye makeup has the eﬀect of making the eyes appear
larger and emphasizing them so that it can improve theface attractive facial impression. Killian et al.7) inves-

170

tigated the change in attractiveness and symmetry, andthen they clariﬁed the correlation between the two by applying makeup to speciﬁc areas of the face. The results
showed that the upper eyelids makeup signiﬁcantly increased attractiveness. Makeup on the areas of increased
attractiveness slightly decreased asymmetry. Thus, theymentioned that makeup may increase attractiveness atleast partially by decreasing asymmetry. Jones et al.8)
veriﬁed the hypothesis that makeup makes facial features(eyebrows, eyes, nose, and mouth) appear larger. The
results showed that both self-applied and professionallyapplied makeup was perceived to have enlarged the eyes.
Morikawa et al.9) revealed through experiments using psychophysical methods that the eyes with eye shadow appeared larger. They showed that the perception of eye
size with eye shadow was related to both the assimilationeﬀect of the eye and the eye shadow, and the assimilation eﬀect of the eyebrow and the eye caused by the eye
shadow ﬁlling the area between the eye and the eyebrow.
These results support that changing the size of facial features with makeup increases facial attractiveness.
Moreover, eyeliner which is a makeup technique drawing a line around the eyelash line has the eﬀect of making
the eyes appear larger. Matsushita et al.10) indicated a
psychological explanation for this eﬀect. They experimentally revealed that the eyeliner had an illusion eﬀectof making the eyes appear larger based on the Delboeuf
illusion, and that the color and the range of eyeliner (“up-IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

per eyelid only” or “upper eyelid + lower eyelid”) aﬀect
how much the eyes appear to be larger (called the amountof illusion). The Delboeuf illusion is the illusion that the
inner circle of a double circle appears to be larger thanthat of a single circle (the inner circle is the same sizeas the single circle), as shown in Fig. 1. This illusion iscaused by the assimilation of the outer and inner circlesof the double circles. Okuri et al.11) considered that the
thickness of eyeliner as well as the color and the rangeof eyeliner which Matsushita et al.10) focused on aﬀects
the amount of illusion, so they measured the amount of
illusion between the thickness of eyeliner and the size ofthe double eyelids using a psychophysical measurement
method. Also, they conducted a questionnaire to investigate the sense of incongruity for the eyes with eyeliner ofvarious thicknesses to determine the appropriate thickness of eyeliner for eye makeup. However, their study
only examined the case of double eyelids, so they have noﬁndings in the case of diﬀerent eyelids shapes.
Eyelids shapes of Asian women are classiﬁed into single, hooded and double, and the attractiveness diﬀers depending on the eyelids shapes. East and Southeast Asian
women often have single eyelids or thin double eyelids
congenitally. Women with double eyelids are sometimes
said to be more attractive than single eyelids because
their eyes appear larger. Jung et al.12) measured thephysical width and height of the eyes in Korean subjects
with and without double eyelids to investigate whether
the eyelids shapes actually make a diﬀerence in eye size.
Their measurement results revealed that the width andheight of the eyes were signiﬁcantly longer in the double
eyelids than in the single eyelids. There were a couple ofreports that the eﬀect of eye makeup diﬀered dependingon the shape of eyelids. Takehara et al.13) examined thephysical and mental attractiveness of woman s faces withsingle/double eyelids and with/without eye bag makeup.
Their results indicated that the faces with double eyelids were more attractive than single eyelids. The faces
with eye bag makeup were generally perceived as less attractive than those without them. Nishino et al.14) examined visual eﬀects of double eyelids with psychologicalexperimental methods and compared with the eﬀects byeyeliner. They showed that the eye enlargement eﬀects
of double eyelids and eyeliner have additive eﬀects. The
experiment suggested that double eyelids are more attractive than single eyelids, and this eﬀect is greater witheyeliner makeup.

171

Fig. 1 An example of the Delboeuf illusion

We have brieﬂy reviewed the studies on attractive
makeup. To our knowledge, there have been no related
studies that investigated the relations the length and thethickness of the attractive eyeliner by eyelids shape (hereafter, referred to as “eyeliner style”). In this paper, wequantitatively measure the length and the thickness ofeyeliner that makes the eyes appear the largest without
any sense of incongruity for each eyelids shape (single eyelids, hooded eyelids, and double eyelids) by continuously
changing the numerical values of both the length and thethickness of the eyeliner. We numerically clarify which
length and thickness of the eyeliner make the eyes look
the largest without any sense of incongruity is an attractive eyeliner style since larger eyes are said to be more
attractive.

2.

MethodWe conducted an experiment to quantitatively verify
the eyeliner style that makes the eyes perceived to be thelargest without any sense of incongruity and considered
whether the eyeliner style has the eﬀect of increasing attractiveness.
2.1

Participants

Thirty-ﬁve people (mean age 22.5 years, SD = 1.79,
22 participants with experienced makeup and 13 participants with inexperienced makeup) with normal visualacuity and color vision participated in the experiment.
All experienced participants in this experiment were female, and all inexperienced participants were male. We
explained the contents of our experiment in verbal andwritten, then participants consented to the experimentof their own free will. Additionally, we used statisticallyprocessed data for analysis so that individuals cannot be
identiﬁed.
2.2

Environment

During the measurement, we instructed that cold lighting should be used, the brightness should be constant,and the direction of the lighting should be ﬁxed. The display was s2430L dell monitor (24 inches). The positionIIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 2 A model of a single eyelids eyeliner

and orientation of the display should be ﬁxed, and theluminance should be measured with a color luminancemeter (CS-100A, Konica Minolta) to be approximately
250 cd/m2. The aspect ratio of the display was 16:9.
The height of the display screen was about 30.5 cm, andthe observation distance was set to about 90 cm (threetimes the height of the display screen). All participants
followed the experimental procedures.
2.3

Eyeliner models

In order to represent eyeliner digitally, we deﬁne eyeliner models based on general shape of eyeliner. The eyeliner models were determined based on the coordinatesof the feature points on the upper eyelids. We use MediaPipe Face Mesh15) to get the feature points. MediaPipe
Face Mesh is a library provided by Google that estimates468 3D facial landmarks in real-time. A model of theeyeliner shape for each eyelids shape is shown as follows.
The models parameterized thickness and length.
Single eyelids have a structure in which the skin is notconnected to the muscles that raise the eyelids, so thatthe skin over the eyelids is not pulled up and covers theeyelids when the eyes are open. Hooded eyelids have asimilar structure on the inner side of the eye. For thisreason, eyeliner models were created for the single eyelids
and the hooded eyelids as far as they are visible when theeyes are open.

172

2.3.1 Single eyelids eyeliner
The model of a single eyelids eyeliner is shown in Fig. 2and it is determined by the following steps.
1. Determine sep1, sep2, sep3, and sep4 from the eye
landmarks．
2. Calculate a sequence of points following the eye
shape by quadratic B-spline interpolation with sep1,
sep2, sep3, and sep4 as control points. Note that thevertical coordinates of the point sequence should be
diﬀerent by 1px.
3. Let sep6 be a point in the point sequence calculatedin Step 2 that matches a vertical coordinate shiftedvertically from sep1 by thickness.
4. The horizontal coordinate of sep5 is the point shiftedhorizontally from the lateral canthus landmark bylength, and the vertical coordinate is the point shiftedvertically from sep1 by thickness (same vertical coordinate as sep6 ).
5. Create a polygon with sep1, sep5, sep6, and interpolation points between sep6 and sep1.
6. Fill in the polygon in Step 5.
2.3.2 Hooded eyelids eyeliner
The model of a hooded eyelids eyeliner is shown in
Fig. 3 and it is determined by the following steps.
1. Determine hep1 hep2, hep3, and hep4 from the eye
landmarks.
2. The point that is shifted length parallel horizontallyfrom the lateral canthus landmark is hep5.
3. hep8 and hep7 are the same points as hep2 and hep3,
respectively.
4. Let hep6 be a point that is shifted vertically from
hep4 by thickness.
5. Create a polygon by quadratic B-spline interpolation
with hep1, hep2, hep3, hep4, hep5, hep6, hep7, andhep8 as control points.
6. Fill in the polygon in Step 5.
2.3.3 Double eyelids eyeliner
The model of a double eyelids eyeliner is shown in
Fig. 4 and it is determined by the following steps.
1. Determine dep1, dep2, dep3, and dep4 from the eye
landmarks.
2. The point that is shifted length parallel horizontallyfrom the lateral canthus landmark is dep5.
3. Let dep8 be a point that is shifted vertically fromdep2 by thickness. Similarly, determine dep7 from
dep3, and dep6 from dep4.
4. Create a polygon by quadratic B-spline interpolation
with dep1, dep2, dep3, dep4, dep5, dep6, dep7, andIIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 3 A model of a hooded eyelids eyeliner

Fig. 4 A model of a double eyelids eyeliner

dep8 as control points.
5. Fill in the polygon in Step 4.
2.4

SystemFigure 5(a) shows the image used in the experiment,a 669 px × 946 px face image of a Japanese woman withdouble eyelids without makeup. The face image was arranged with hair style that covered the eyebrows, because
the eyebrows and hair may cause bias in the size perception. The right eye was created by inverting the left eyethen the eyes were placed symmetrically in a sagittal line.
We retouched the boundary of the right eye and the facial skin by Photoshop to make them natural. We created
two other face images with single eyelids and hooded eyelids by retouching the face image with double eyelids asshown in Fig. 5(b)(c). Since the eye size of single eyelids
is physically shorter than double eyelids in both widthand height12), the face image with single eyelids was created to have smaller eyes than the face image with double
eyelids. The hooded eyelids were created by erasing the
double eyelids line from the inner corner to the center ofthe eye without changing the size of the double eyelids. In
this way, three models, a double eyelids model, a hooded
eyelids model, and a single eyelids model were prepared.
The actual image on the screen was set to 1 px = 0.3
mm (rounded to the second decimal place) at a resolutionof 1,920 px × 1,080 px. The eye width on the screen
was adjusted to be approximately 2.7 cm and the eyeheight was approximately 1.0 cm in the hooded eyelids
model and the double eyelids model based on the averageinterocular breadth of 3.7 cm and biectocanthion breadthof 9.0 cm, which are typical for Japanese young adultfemale16). For the single eyelids model, the width was173

approximately 2.5 cm and the height was approximately
0.8 cm.
In this experiment, we prepared a system that canchange the parameters of length and thickness of the eyeliner drawn by computer graphics simulation to quantitatively verify the eyeliner style. The system allows theparticipants to adjust length and thickness by using theslider interface. The adjustable range is 0 mm to about
11.1 mm (0 px to 40 px) for length of three types of eyelids shapes models, 0 mm to about 4.2 mm (0 px to 15px) for thickness of the single eyelids model, and 0 mm toabout 5.5 mm (0 px to 20 px) for thickness of the hooded
eyelids model and the double eyelids model. The eyeliner was drawn for each eyelids shape using the eyeliner
models deﬁned in Section 2.3.
Some of the eyeliner for each eyelids shape model isshown in Table 1 for length manipulation and in Table 2
for thickness manipulation. In these ﬁgures, only the left
eye is cropped. The images displayed in the experimentwere the entire face.
2.5

Procedure

Participants looked at the images presented by the experimental system and adjusted length and thickness ofthe eyeliner to get the most appropriate eyeliner for eachcondition.
The experimental procedure was as follows.
First, the screen showed an image of a face whose thelength and the thickness of eyeliner can be changed sequentially by the participant’s manipulation.
Then, the participants were asked to respond to three
conditions per model (the single eyelids model, thehooded eyelids model, and the double eyelids model) byIIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）(a) Double eyelids
Fig. 5(b) Single eyelids(c) Hooded eyelids

The face images used in the experiment with diﬀerent eyelids shapes
Table 1 Eyeliner for length manipulation

Table 2 Eyeliner for thickness manipulation

adjusting only length or only thickness for a total of six
requirements. The following three conditions were given:
1. make the eyeliner that makes the eyes appear thelargest, 2. make the eyeliner at the limit of the senseof incongruity, and 3. make the face the most attractivewith eyeliner.
We did not set time limits in the experiment, and wesuggested that the participants take breaks at their own
discretion so that they could make correct decisions. We
also instructed them to answer the questions as intuitively
as possible.

174

The trial was conducted in the order of the single eyelids model, the hooded eyelids model, and the double eyelids model. We conducted these trials for three models
and provided a short break（at least one-minute） beforechanging models and conducting the next trial.
After conducting the veriﬁcation experiment, we surveyed the participants about the criteria for the sense ofincongruity and attractiveness. The survey asked for text
open-ended questions such as “In which cases did you feelthe sense of incongruity?” and “In which cases did you
ﬁnd most attractive?”IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

In this experiment, we ﬁxed the thickness for adjusting only the length of the eyeliner as the standard thickness and the length for adjusting only the thickness as
the standard length. The standard thickness was ﬁxed
at 2.8 mm for the single eyelids model, 1.4 mm for the
hooded eyelids model and the double eyelids model. For
the hooded eyelids model and the double eyelids model,
the standard thickness was determined based on 1.5 mm
brush tip of a typical ultraﬁne pencil eyeliner. For the
single eyelids model, the standard thickness was determined based on the results of a preliminary experimentconducted by three participants with makeup experience
who were not part of the experiment because the shapeof the eyeliner is quite diﬀerent from other models. The
standard length was ﬁxed at 2.8 mm for three models.
For the double eyelids model, we determined the standard length so that the horizontal coordinate was equal
to the point at the end of the double eyelid line on theouter corner of the eye. We adopted the same standard
length for the single eyelids model and the hooded eyelids
model.

3.
3.1

Results and Analysis

Results

The experimental results for the length and the thickness of eyeliner are shown in Fig. 6 to Fig. 8 by each eyelids shape using the box plot. The vertical axis indicatesthe length or thickness of the eyeliner, and the horizontal
axis indicates the condition. The larger the value of eyeliner length or thickness, the heavier the eye makeup, andthe smaller the value, the lighter the eye makeup. Cross
marks in the graph indicate mean values. Then, the meanvalues for each condition are shown in Table 3 by eacheyelids shape. Moreover, Fig. 9 to Fig. 11 show theresults of the experiment separately for the experienced
makeup and the inexperienced makeup participants. In
Table 3 and Fig. 6 to Fig. 11, “Largest” indicates theeyeliner that makes the eyes appear the largest, “Incongruity” indicates the eyeliner at the limit of the sense ofincongruity, and “Attractive” indicates the eyeliner withthe most attractive. In Table 3, “Largest & Incongruity”
shows the smaller value between “Largest” and “Incongruity.”
3.2

Analysis

The purpose of this experiment is to verify whether
the eyeliner style that makes the eyes appear the largest
without any sense of incongruity is an attractive eyeliner

175

style. The eﬀect of the eye size illusion becomes larger
when the thickness of the eyeliner is thicker. There havebeen found that larger eyes make attractive, but excessively thick eyeliner is not appropriate as eye makeup
because it has any sense of incongruity11). Therefore,
we conducted an equivalence test to conﬁrm whether theeyeliner style that makes the eyes appear the largest without any sense of incongruity and the eyeliner style withthe most attractive are equivalent. An equivalence test
demonstrates statistical equivalence. Equivalence can beproven if the 95% conﬁdence interval of the diﬀerence between levels is within a predeﬁned equivalence margin (arange in which no signiﬁcant diﬀerence is determined).
In this experiment, the equivalence margin was set to 1.5mm for both length and thickness because the brush tipof a typical ultra-ﬁne pencil eyeliner is 1.5 mm17). Figure 12 shows the results of the equivalence test with theequivalence margin and 95% conﬁdence intervals for thediﬀerence between the means of “Largest & Incongruity”
and “Attractive” for each eyelids shape. The 95% conﬁdence intervals for all eyelids shapes for both eyeliner
length and thickness were within the equivalence margin.
The results of the equivalence test showed that there wasno signiﬁcant diﬀerence between “the eyeliner style thatmakes the eyes appear the largest without any sense ofincongruity” and “the eyeliner style with the most attractive” for each eyelids shape and indicated that they
were equivalent. As a combination of the standard length
or the standard thickness, the length of 6.3 mm and thethickness of 2.6 mm for single eyelids, the length of 6.6mm and the thickness of 3.4 mm for hooded eyelids, andthe length of 7.0 mm and the thickness of 2.6 mm for
double eyelids were shown to make the eyes appear thelargest without any sense of incongruity and also to makethe face more attractive (Table 4).
The equivalence test was also conducted on both theexperienced and the inexperienced participants. The results showed that “the eyeliner style that makes the eyes
appear the largest without any sense of incongruity” and“the eyeliner style with the most attractive” except for
the length of double eyelids on the experienced participants were equivalent. For the length of double eyelids
on the experienced participants, 95% conﬁdence interval
was out of the equivalence margin (mean diﬀerence 0.74,
95%CI [-0.11, 1.59]).

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 6 A result of the single eyelids eyeliner

Fig. 7 A result of the hooded eyelids eyeliner

Fig. 8 A result of the double eyelids eyeliner

Table 3 A result of the experiment (mean value)176

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

Fig. 9 A result of the single eyelids eyeliner (Comparison of experienced makeup
and inexperienced makeup participants)

Fig. 10 A result of the hooded eyelids eyeliner (Comparison of experienced makeup
and inexperienced makeup participants)

Fig. 11 A result of the double eyelids eyeliner (Comparison of experienced makeup
and inexperienced makeup participants)

4.

Discussion

4.1 Attractive eyeliner
The equivalence test of the experimental results showedthat the eyeliner which makes the eyes appear the largest

177

without any sense of incongruity is equivalent to the most
attractive eyeliner. These results suggest that the eyeliner
style that can make the eyes appear the largest within notexcessive as makeup is attractive. The post experimentIIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）
Table 4 Eyeliner that made the eyes appear the largest
without any sense of incongruity

Fig. 12 A result of the equivalence test

survey on attractiveness criteria supported the results ofthis veriﬁcation experiment since more than half of the respondents felt that the faces with moderately large eyes
were attractive. However, the length of eyeliner whichmakes the eyes appear the largest without any sense of incongruity on the double eyelids was not equivalent to themost attractive length of eyeliner for the experienced participants. In the case of the double eyelids, the length of
the attractive eyeliner seemed to reﬂect preference rather
than making the eyes appear larger.
The eyes with excessive eyeliner make the boundary
between the face and the eyes clear, and the geometric
illusion caused by the contrast between the surrounding
the face and the eyes conversely makes the eyes appear
smaller10). According to the experiment of measuring theeye size illusion by the thickness of eyeliner11), the eyes
appeared to be larger up to a certain level of thickness ofthe eyeliner, but the illusion eﬀect decreased beyond thatlevel. In our experiment, the peak of attractiveness andthe peak of appearing to have larger eyes without any
sense of incongruity were equivalent. These results suggest that excessive eyeliner not only decreases the illusionof larger eyes but also may not be attractive.
4.2

Diﬀerences for eyelids shape

Since the eyeliner model diﬀers signiﬁcantly between
the single and the hooded/double eyelids, we discuss
the results separately for the single eyelids and thehooded/double eyelids.
4.2.1

Single eyelids

Figure 9 shows that the IQR are larger for the inexperienced participants than for the experienced participants,
indicating that there are larger scatter in the data. The

178

experienced participants are familiar with eyeliner as a
makeup method, and it is easy for them to set a standard
based on the eyeliner they usually wear. However, theinexperienced participants are not familiar with eyeliner,
and it is diﬃcult for them to set a standard. In particular, for the single eyelids eyeliner, the eyeliner style wasmore unfamiliar than other eyeliner styles because it was
drawn only on the outer corners of the eyes. This is a reason for the large scatter in the data for the inexperienced
participants.
4.2.2

Hooded eyelids and double eyelids

Comparing the results for hooded eyelids with double eyelids (Fig. 7 and Fig. 8), there was a signiﬁcantdiﬀerence in “Attractiveness” for the thickness of eyeliner between the double eyelids and the hooded eyelids(Wilcoxon signed rank test, p <.01, 95%CI[0.69, 1.11]),
with smaller values for the double eyelids. This suggeststhat the light eyeliner was preferred for double eyelids.
In terms of an illusion, the double eyelids alone have an
illusion eﬀect of making the eyes appear larger by around
4.5%10),14). Furthermore, the combination of eyeliner andthe double eyelids produces an illusion eﬀect of nearly
6%14). Since the eﬀect of the illusion on the face9)–11) is
not so great compared to a geometric illusion of simple
circle combinations, it can be assumed that the eﬀect ofthe makeup illusion on the face is limited to the extent
that it does not become unnatural, and no more the eﬀect
can be obtained. The double eyelids already have clear
double eyelids line in the vertical direction, so even if the
thickness of eyeliner is more modest than the hooded eyelids, the illusion eﬀect could be enough to make the facelook more attractive. The post experiment survey showed

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）

that more than half of the respondents felt the sense ofincongruity when the eyeliner was thicker than the double eyelids line for hooded and double eyelids. Moreover,some respondents felt that even light eyeliner (eyeliner ofmodest length and thickness) that clearly shows the double eyelids line for double eyelids was attractive. Therefore, the double eyelids line is an important factor to determine the eyeliner style.
The thickness of eyeliner that makes the eyes appear
the largest without any sense of incongruity for double
eyelids in this experiment was 2.6 mm. According to
the psychophysical measurement11), it was about 2.4 mm.
Since the diﬀerence of 0.2 mm is less than 1 px, the resultsof this experiment is approximated as the psychophysical measurement result. However, while Okuri et al.11)
reported diﬀerences in the illusion eﬀect between maleand female participants depending on the thickness ofeyeliner, there was no diﬀerence between male (the inexperienced participants) and female (the experienced participants) in the “Largest” for the thickness (Fig. 11) ofeyeliner in our experiment (Wilcoxon rank sum test, p= .14, 95%CI[-1.66, 0.27]) . In their study, the illusioneﬀect increased with thicker eyeliner only for male, while
the illusion eﬀect did not vary with thicker eyeliner for
female. In contrast, female also perceived that the thickeyeliner made the eyes appear larger in our experiment,
and the eyeliner that makes the eyes appear the largest
was approximately the same for both male and female.
We assume that this diﬀerence in results is due to thediﬀerence in experimental methods. Their experimentused the staircase method to discretely present diﬀerent
level of the length and thickness of the eyeliner. In our approach for using the slider, the user interactively and continuously changes the stimulus so the approach can measure more detailed makeup diﬀerences. Moreover, in theirexperiment, the image presentation time was as short as
3 seconds, whereas in our experiment, the presentation
time was not limited, and the participants adjusted theslider while looking at the image carefully. These diﬀerences in experimental methods may have some aﬀected
the results.

5.

Conclusions

In this paper, we conducted a veriﬁcation experimentin which we quantitatively investigated the eyeliner stylethat can make the eyes appear the largest without any
sense of incongruity and the most attractive eyeliner stylefor each eyelids shape by using the eyeliner model. The

179

results of the experiment indicated that the eyeliner thatmakes the eyes appear the largest without any sense ofincongruity is equivalent to the most attractive eyeliner.
Furthermore, excessive eyeliner was found to decrease attractiveness. These ﬁndings are a contribution to thestudy of digital makeup.
The eyeliner models were determined from general eyeliner shapes. Therefore, the result obtained in this experiment is an attractive eyeliner in the case of the generalshape model. Further experiments are needed in orderto obtain ﬁndings on other eyeliner shapes. The combination of the length and the thickness in cases where
no standard length or thickness has been determined is a
future issue. Consideration of eye characteristics such aseye aspect ratio and iris ratio is a future topic.
Eye makeup includes not only eyeliner but also eyeshadow and mascara. Generally, if people wear heavy
makeup with thick eyeliner, they also wear other makeup
as heavy. Therefore, it is important to consider the total balance of the whole eye makeup. It is a future task
to develop methods that can adjust other eye makeup,
such as eye shadow and mascara, with parameters and toinvestigate the total balance of the eye makeup.
In the future, we will consider developing a system thatcan recommend eye makeup styles that increase attractiveness according to the individual’s eyelids shape basedon the results obtained in this paper. Moreover, obtaining feedback from makeup experts would be future tasks.

References
1) R. L. Dipboye, R. D. Arvey, D. E. Terpstra: “Sex and PhysicalAttractiveness of Raters and Applicants as Determinants ofResumé Evaluations”, Journal of Applied Psychology, Vol.62,
No.3, pp.288-294 (1977).
2) J. Y. Baudouin, G. Tiberghien: “Symmetry, Averageness, andFeature Size in the Facial Attractiveness of Women”, Acta Psychologica, Vol.117, No.3, pp.313-332 (2004).
3) D. Jones, C. L. Brace, W. Jankowiak, K. N. Laland, L. E.
Musselman, J. H. Langlois, L. A. Roggman, D. Pérusse, B.
Schweder, D. Symons: “Sexual Selection, Physical Attractiveness and Facial Neoteny: Cross-Cultural Evidence and Implications.”, Current Anthropology, Vol.36, No.5, pp.723-748(1995).
4) D. I. Perrett, K. J. Lee, Penton-Voak, D. Rowland, S.
Yoshikawa, D. M. Burt, S. P. Henzi, D. L. Castles, S. Akamatsu: “Eﬀects of Sexual Dimorphism on Facial Attractiveness”, Nature, Vol.394, No.6696, pp.884-887 (1998).
5) K. Tagai, H. Ohtaka, H. Nittono: “Faces with Light Makeup
Are Better Recognized than Faces with Heavy Makeup”, Frontiers in Psychology, Vol.7, No.226 (2016).
6) R. Mulhern, G. Fieldman, T. Hussey, J. Lévêque, P. Pineau:“Do Cosmetics Enhance Female Caucasian Facial Attractiveness?”, International Journal of Cosmetic Science, Vol.25,
No.4, pp.199-205 (2003).

IIEEJ Transactions on Image Electronics and Visual Computing Vol.10 No.2 （2022）
7) A. C. Killian, S. Mitra, J. J. Peissig: “The Role of Regional
Contrast Changes and Asymmetry in Facial Attractiveness
Related to Cosmetic Use”, Frontieres in Psychology, Vol.9,
No.2448, pp.1-11 (2018).
8) A. L. Jones, A. Porcheron, R. Russell: “Makeup Changes the
Apparent Size of Facial Features”, Psychology of Aesthetics,
Creativity, and the Arts, Vol.12, No.3, pp.359-368 (2018).
9) K. Morikawa, S. Matsushita, A. Tomita, H. Yamanami: “A
Real-Life Illusion of Assimilation in the Human Face: Eye Size
Illusion Caused by Eyebrows and Eye Shadow”, Frontiers in
Human Neuroscience, Vol.9, No.139, pp.1-9 (2015).
10) S. Matsushita, K. Morikawa, H. Yamanami: “Measurementof Eye Size Illusion Caused by Eyeliner, Mascara, and Eye
Shadow”, Journal of Cosmetic Science, Vol.66, No.3, pp.161174 (2015).
11) R. Okuri, S. Kodama, T. Takahashi: “Measurement of Eye
Size Illusion Caused by Thickness of Eyeliner on Double-Eyelid
Eyes”, IIEEJ Transactions on Image Electronics and VisualComputing, Vol.10, No.1, pp.47-54 (2022).
12) H. B. Jung, D. G. Han, J. S. Shim, Y. J. Lee, S. Kim: “Comparison of Eye Measurements between Young Korean Women
with Inborn Double Eyelids and Those with Single Eyelids”,
Archives of Aesthetic Plastic Surgery, Vol.26, No.1, pp.7-11(2020).
13) T. Takehara, N. Inoue, L. Yamamoto, M. Shimizu: “Faces with
Double Upper Eyelids are Attractive while Lower Eyelids with
Eye Bag Makeup might be Unattractive”, Journal of Japan
Society of Kansei Engineering, Vol.20, No.2, pp.121-128 (2021)(in Japanese).
14) K. Nishino, N. Okiyama: “A Psychological Study for the Visual Eﬀects Induced by Double Eyelids - Compared with Eyeliner and Considered for Gender Diﬀerence”, FRAGRANCEJOURNAL, Vol.49, No.5, pp.46-52 (2021) (in Japanese).
15) C. Lugaresi, J. Tang, H. Nash, C. McClanahan, E. Uboweja, M.
Hays, F. Zhang, C. Chang, M. G. Young, J. Lee, W. Chang, W.
Hua, M. George, M. Grundmann: “MediaPipe: A Frameworkfor Building Perception Pipelines”, arXiv (2019).
16) M. Kouchi, M. Mochimaru: “Anthropometric Database ofJapanese Head 2001”, National Institute of Advanced Industrial Science and Technology, H16PRO-212 (2008).
17) CANMAKE TOKYO, Creamy Touch Liner, https://www.
canmake.com/en/item/detail/36 (2022).

(Received September 9, 2022)(Revised October 29, 2022)

180

Rena Okuri
She received her B.E. and M.E degree from
Tokyo Denki University in 2020 and 2022.
Her research interests include facial image
processing.

Shuhei Kodama
He received his Ph.D. (Engineering) from
Tokyo Denki University in 2019. He was anengineer at ASTRODESIGN, Inc. from 2019
to 2021. He is currently an assistant professor at Tokyo Denki University. His researchﬁelds are NPR, HCI, and fabrication.

Tokiichiro Takahashi (Fellow)
He received his B.E. degrees from Niigata
Uiversity in 1977, and Ph. D. from the University of Tokyo in 2005, respectively. He
had been working at Nippon Telegraph andTelephone Corporation and engaged in Pattern Recognition, Computer Graphics, andLearning Science since 1977. Since 2003, hehas been a professor at Tokyo Denki University. He has been a visiting researcher at ASTRODESIGN,Inc. since 2017. His researchinterests include computer graphics, image
processing, and visual computing.

Revised: January 6, 2017
Revised: July 6, 2018

Guidance for Paper Submission
1.

Submission of Papers(1) Preparation before submission
・ The authors should download “Guidance for Paper Submission” and “Style Format” from the“Academic Journals”, “English Journals” section of the Society website and prepare the paperfor submission.
・ Two versions of “Style Format” are available, TeX and MS Word. To reduce publishing costs
and effort, use of TeX version is recommended.
・ There are four categories of manuscripts as follows: Ordinary paper: It should be a scholarly thesis on a unique study, development orinvestigation concerning image electronics engineering. This is an ordinary paper topropose new ideas and will be evaluated for novelty, utility, reliability andcomprehensibility. As a general rule, the authors are requested to summarize a paper
within eight pages.
 Short paper: It is not yet a completed full paper, but instead a quick report of the partial
result obtained at the preliminary stage as well as the knowledge obtained from the said
result. As a general rule, the authors are requested to summarize a paper within four
pages.
 System development paper: It is a paper that is a combination of existing technology or it
has its own novelty in addition to the novelty and utility of an ordinary paper, and thedevelopment results are superior to conventional methods or can be applied to othersystems and demonstrates new knowledge. As a general rule, the authors are requested
to summarize a paper within eight pages.
 Data Paper: A summary of data obtained in the process of a survey, product development,
test, application, and so on, which are the beneficial information for readers even though
its novelty is not high. As a general rule, the authors are requested to summarize a paper
within eight pages.
・ To submit the manuscript for ordinsry paper, short paper, system development paper, or data
paper, at least one of the authors must be a member or a student member of the society.
・ We prohibit the duplicate submission of a paper. If a full paper, short paper, systemdevelopment paper, or data paper with the same content has been published or submitted toother open publishing forums by the same author, or at least one of the co-authors, it shallnot be accepted as a rule. Open publishing forum implies internal or external books,
magazines, bulletins and newsletters from government offices, schools, company
organizations, etc. This regulation does not apply to a preliminary draft to be used at anannual meeting, seminar, symposium, conference, and lecture meeting of our society or othersocieties (including overseas societies). A paper that was once approved as a short paper andbeing submitted again as the full paper after completion is not regarded as a duplicate
submission.
(2) Submission stage of a paper・ Delete all author information at the time of submission. However, deletion of reference
information is the author’s discretion.
・ At first, please register your name on the paper submission page of the following URL, andthen log in again and fill in the necessary information. Use the “Style Format” to upload yourmanuscript. An applicant should use PDF format (converted from dvi of TeX or MS Word

181

format) for the manuscript. As a rule, charts (figures and tables) shall be inserted into themanuscript to use the “Style Format”. (a different type of data file, such as audio and video,
can be uploaded at the same time for reference.)
http://www.editorialmanager.com/iieej/
・ If you have any questions regarding the submission, please consult the editor at our office.
Contact:
Person in charge of editing
The Institute of Image Electronics Engineers of Japan
3-35-4-101, Arakawa, Arakawa-Ku, Tokyo 116-0002, Japan
E-mail: hensyu@iieej.org
Tel: +81-3-5615-2893, Fax: +81-3-5615-2894

2.

Review of Papers and Procedures(1) Review of a paper・ A manuscript is reviewed by professional reviewers of the relevant field. The reviewer willdeem the paper “acceptance”, “conditionally acceptance” or “returned”. The applicant is
notified of the result of the review by E-mail.
・ Evaluation method
Ordinary papers are usually evaluated on the following criteria:
 Novelty: The contents of the paper are novel.
 Utility: The contents are useful for academic and industrial development.
 Reliability: The contents are considered trustworthy by the reviewer.
 Comprehensibility: The contents of the paper are clearly described and understood by
the reviewer without misunderstanding.
Apart from the novelty and utility of an ordinary paper, a short paper can be evaluated byhaving a quickness on the research content and evaluated to have new knowledge withresults even if that is partial or for specific use.
System development papers are evaluated based on the following criteria, apart from thenovelty and utility of an ordinary paper.
 Novelty of system development: Even when integrated with existing technologies, thenovelty of the combination, novelty of the system, novelty of knowledge obtained fromthe developed system, etc. are recognized as the novelty of the system.
 Utility of system development: It is comprehensively or partially superior compared tosimilar systems. Demonstrates a pioneering new application concept as a system. The
combination has appropriate optimality for practical use. Demonstrates performance
limitations and examples of performance of the system when put to practical use.
Apart from the novelty and utility of an ordinary paper, a data paper is considered novel if
new deliverables of test, application and manufacturing, the introduction of new technology
and proposals in the worksite have any priority, even though they are not necessarily
original. Also, if the new deliverables are superior compared to the existing technology andare useful for academic and industrial development, they should be evaluated.
(2) Procedure after a review
・ In case of acceptance, the author prepares a final manuscript (as mentioned in 3.).
・ In the case of acceptance with comments by the reviewer, the author may revise the paper inconsideration of the reviewer’s opinion and proceed to prepare the final manuscript (as182

mentioned in 3.).
・ In case of conditional acceptance, the author shall modify a paper based on the reviewer’s
requirements by a specified date (within 60 days), and submit the modified paper for
approval. The corrected parts must be colored or underlined. A reply letter must be attachedthat carefully explains the corrections, assertions and future issues, etc., for all of theacceptance conditions.
・ In case a paper is returned, the author cannot proceed to the next step. Please look at thereasons the reviewer lists for the return. We expect an applicant to try again after reviewing
the content of the paper.
(3) Review request for a revised manuscript
・ If you want to submit your paper after conditional acceptance, please submit the reply letterto the comments of the reviewers, and the revised manuscript with revision history to thesubmission site. Please note the designated date for submission. Revised manuscripts
delayed more than the designated date be treated as new applications.
・ In principle, a revised manuscript willl be reviewed by the same reviewer. It is judged either
accceptance or returned.
・ After the judgment, please follow the same procedure as (2).

3.

Submission of final manuscript for publication(1) Submission of a final manuscript
・ An author, who has received the notice of “Acceptance”, will receive an email regarding thecreation of the final manuscript. The author shall prepare a complete set of the final
manuscript (electronic data) following the instructions given and send it to the office by thedesignated date.
・ The final manuscript shall contain a source file (TeX edition or MS Word version) and a PDF
file, eps files for all drawings (including bmp, jpg, png), an eps file for author’s photograph(eps or jpg file of more than 300 dpi with length and breadth ratio 3:2, upper part of the body)for authors’ introduction. Please submit these in a compressed format, such as a zip file.
・ In the final manuscript, write the name of the authors, name of an organizations,introduction of authors, and if necessary, an appreciation acknowledgment. (cancel macros in
the Style file)・ An author whose paper is accepted shall pay a page charge before publishing. It is theauthor’s decision to purchase offprints. (ref. page charge and offprint price information)(2) Galley print proof
・ The author is requested to check the galley (hard copy) a couple of weeks before the paper ispublished in the journal. Please check the galley by the designated date (within one week).
After making any corrections, scan the data and prepare a PDF file, and send it to our officeby email. At that time, fill in the Offprint Purchase Slip and Copyright Form and return thescanned data to our office in PDF file form.
・ In principle, the copyrights of all articles published in our journal, including electronic form,
belong to our society.
・ You can download the Offprint Purchase Slip and the Copyright Form from the journal on our
homepage. (ref. Attachment 2: Offprint Purchase Slip, Attachment 3: Copyright Form)(3) Publication
・ After final proofreading, a paper is published in the Academic journal or English transaction(both in electronic format) and will also be posted on our homepage.

183

IIEEJ Trans. on Image Electronics and Visual Computing©2022 IIEEJ

Vol. 10, No. 2 2022

Editor in Chief: Osamu Uchida
The Institute of Image Electronics Engineers of Japan
3-35-4-101，Arakawa，Arakawa-ku，Tokyo 116-0002，Japan