IIEEJ Transactions on
IIEEJ Transactions on
Image Electronics and Visual Computing

Image Electronics andVisual Computing

Special Issue on Journal Track Papers in IEVC2021

Vol. 9, No. 2 2021

I.I.E.E.J

Editorial Committee of IIEEJ                              (GLWRULQ&KLHI0HL.2'$0$ +LURVKLPD8QLYHUVLW\ 

9LFH(GLWRUVLQ&KLHI
2VDPX8&+,'$ 7RNDL8QLYHUVLW\  
1DRNL.2%$<$6+, 6DLWDPD0HGLFDO8QLYHUVLW\ <XULNR7$.(6+,0$ 7RN\R8QLYHUVLW\RI7HFKQRORJ\ $GYLVRU\%RDUG<DVXKLNR<$68'$ :DVHGD8QLYHUVLW\(PHULWXV  +LGH\RVKL720,1$*$ :DVHGD8QLYHUVLW\(PHULWXV  .D]XPL.20,<$ .DQDJDZD,QVWLWXWHRI7HFKQRORJ\  )XPLWDND212 7RN\R3RO\WHFKQLF8QLYHUVLW\(PHULWXV  <RVKLQRUL+$725, 7RN\R,QVWLWXWHRI7HFKQRORJ\  
0LWVXML0$7680272 :DVHGD8QLYHUVLW\(PHULWXV  .L\RVKL7$1$.$ 6KLQVKX8QLYHUVLW\ 
6KLJHR.$72 8WVXQRPL\D8QLYHUVLW\(PHULWXV  (GLWRUV<RVKLQRUL$5$, 7RN\R3RO\WHFKQLF8QLYHUVLW\ &KHH6HQJ&+$1 8QLYHUVLW\RI0DOD\D  
1DLZDOD3&+$1'5$6,5, .RJDNXLQ8QLYHUVLW\  &KLQWKDND35(0$&+$1'5$ 6KLEDXUD,QVWLWXWHRI7HFKQRORJ\ 
0DNRWR)8-,6$:$ 8QLYHUVLW\RI7VXNXED ,VVHL)8-,6+,52 .HLR8QLYHUVLW\  .D]XKLNR+$0$0272 7RNDL8QLYHUVLW\  
0DGRND+$6(*$:$ 8WVXQRPL\D8QLYHUVLW\  
5\RVXNH+,*$6+,.$7$ )8-,),/0%XVLQHVV,QQRYDWLRQ&RUS <XNL,*$5$6+, 0HLML8QLYHUVLW\  
7RPRND]X,6+,.$:$㸦7R\R8QLYHUVLW\㸧
0DVDKLUR,6+,.$:$ 㸦6DLWDPD0HGLFDO8QLYHUVLW\㸧
1DRWR.$:$085$ &DQRQ2%  
6KXQLFKL.,085$ )8-,),/0%XVLQHVV,QQRYDWLRQ&RUS 
6KRML.85$.$.( 177'2&202  .D]XWR.$0,.85$ 7RN\R3RO\WHFKQLF8QLYHUVLW\ 
7DNDVKL.$1$, 7KH8QLYHUVLW\RI7RN\R  
7HWVXUR.8*( 1+.(QJLQHHULQJ6\VWHP,QF  .RML0$.,7$ &DQRQ,QF 
7RPRDNL025,<$ 7RN\R'HQNL8QLYHUVLW\  
3DUDPHVUDQ5$9((1'5$1 8QLYHUVLW\RI0DOD\D  .DLVHL6$.85$, ':$1*2&R/WG .RNL6$72 6KRQDQ,QVWLWXWHRI7HFKQRORJ\  
6\XKHL6$72㸦8QLYHUVLW\RI7R\DPD㸧
0DVDQRUL6(.,12 )8-,),/0%XVLQHVV,QQRYDWLRQ&RUS .D]XPD6+,12'$ 8WVXQRPL\D8QLYHUVLW\  
0LNLR6+,1<$ 7RKR8QLYHUVLW\  
6KLQLFKL6+,5$.$:$ $R\DPD*DNXLQ8QLYHUVLW\  .HQLFKL7$1$.$ 1DJDVDNL,QVWLWXWHRI$SSOLHG6FLHQFH  <XNLKLUR768%26+,7$ )XML;HUR[&R/WG 'DLVXNH768'$ 6KLQVKX8QLYHUVLW\ 
0DVDKLUR72<285$ 8QLYHUVLW\RI<DPDQDVKL .D]XWDNH8(+,5$ .DQDJDZD,QVWLWXWHRI7HFKQRORJ\ <XLFKLUR<$0$'$ *HQHVLV&RPPHUFH&R/WG  
1RULPDVD<26+,'$ 1LKRQ8QLYHUVLW\  
7RVKLKLNR:$.$+$5$ )XNXRND,QVWLWXWHRI7HFKQRORJ\2%  .RN6KHLN:21* 0RQDVK8QLYHUVLW\0DOD\VLD  

5HYLHZHU+HUQDQ$*8,55( 6KLQVKX8QLYHUVLW\  .HQLFKL$5$.$:$ 177$GYDQFHG7HFKQRORJ\&RUSRUDWLRQ  
6KRLFKL$5$., 3DQDVRQLF&RUSRUDWLRQ  
7RPRKLNR$5,.$:$ 177(OHFWURQLFV&RUSRUDWLRQ  <XH%$2 7RN\R&LW\8QLYHUVLW\  
1RUGLQ%,15$0/, 0,026%HUKDG  <RRQJ&KRRQ&+$1* 0XOWLPHGLD8QLYHUVLW\  
5RELQ%LQJ<X&+(1 1DWLRQDO7DLZDQ8QLYHUVLW\ .L\RQDUL)8.8( 7RNDL8QLYHUVLW\ 
0RFKDPDG +$5,$', 6HSXOXK1RSHPEHU,QVWLWXWHRI7HFKQRORJ\ 
0DVDNL+$<$6+, 8336$/$8QLYHUVLW\ 
7DNDKLUR+21*8 1(&(QJLQHHULQJ/WG  <XXNRX+25,7$ 8QLYHUVLW\RI7R\DPD  
7DND\XNL,72 2FKDQRPL]X8QLYHUVLW\  
0DVDKLUR,:$+$6+, 1DJDRND8QLYHUVLW\RI7HFKQRORJ\  
0XQHWRVKL,:$.,5, 1DWLRQDO'HIHQVH$FDGHP\RI-DSDQ  <RVKLKLUR.$1$025, 8QLYHUVLW\RI7VXNXED  
6KXQLFKL.$1(.2 +RNNDLGR8QLYHUVLW\ <RXVXQ.$1* 7RN\R3RO\WHFKQLF8QLYHUVLW\ 
3L]]DQX.$121*&+$,<26 &KXODORQJNRUQ8QLYHUVLW\  +LGHWRVKL.$7680$ 7DPD$UW8QLYHUVLW\2%   
0DVDNL.,7$*2 &DQRQ,QF  $NL\XNL.2'$7( 7VXGD&ROOHJH  +LGHNL.20$*$7$ 6DLWDPD0HGLFDO8QLYHUVLW\  <XVKL.20$&+, .RNXVKLNDQ8QLYHUVLW\  
7RVKLKLUR.200$ 7RN\R0HWURSROLWDQ8QLYHUVLW\  
7VXQH\D.85,+$5$ +LWDFKL/WG  
7RVKLKDUX.8526$:$ 0DWVXVKLWD(OHFWULF,QGXVWULDO&R/WG2%  .D]XIXPL.$1('$ +LURVKLPD8QLYHUVLW\  ,WDUX.$1(.2 7RN\R3RO\WHFKQLF8QLYHUVLW\  
7HFN&KDZ/,1* 8QLYHUVLW\RI0DOD\D &KX.LRQJ/22 8QLYHUVLW\RI0DOD\D )
;LDR\DQJ0$2 8QLYHUVLW\RI<DPDQDVKL  .RLFKL0$768'$ ,ZDWH3UHIHFWXUDO8QLYHUVLW\  
0DNRWR0$768., 1774XDULV&RUSRUDWLRQ2%  
7DNHVKL0,7$ 7RVKLED&RUSRUDWLRQ  +LGHNL0,7680,1( 1+.6FLHQFH 7HFKQRORJ\5HVHDUFK/DERUDWRULHV  
6KLJHR025,6+,0$ :DVHGD8QLYHUVLW\  .RXLFKL0876885$ 6KLQV\X8QLYHUVLW\  <DVXKLUR1$.$085$ 1DWLRQDO'HIHQVH$FDGHP\RI-DSDQ  .D]XKLUR12720, .DQDJDZD,QVWLWXWHRI7HFKQRORJ\  
7DNDR212<( 2VDND8QLYHUVLW\  +LGHIXPL26$:$ &DQRQ,QF  .HDW.HRQJ3+$1* 8QLYHUVLW\RI0DOD\D  )XPLKLNR6$,72 *LIX8QLYHUVLW\  
7DNDIXPL6$,72 7RN\R8QLYHUVLW\RI$JULFXOWXUHDQG7HFKQRORJ\ 
7VX\RVKL6$,72 7RN\R,QVWLWXWHRI7HFKQRORJ\  
0DFKLNR6$72 7RN\R3RO\WHFKQLF8QLYHUVLW\(PHULWXV  
7DND\RVKL6(0$6$ 0LWVXELVKL(OHFWULF&RUS2% .DRUX6(=$., 7KH8QLYHUVLW\RI7RN\R  -XQ6+,0$085$ 177  
7RPR\RVKL6+,02%$%$ &KLED8QLYHUVLW\  .DWVX\XNL6+,12+$5$ .RJDNXLQ8QLYHUVLW\  .HLLFKLUR6+,5$, 6KLQVKX8QLYHUVLW\ (LML68*,6$., 1'HVLJQ,QF -DSDQ 'DZQ3XUSOH,QF 3KLOLSSLQHV .XQLKLNR7$.$12 7RN\R0HWURSROLWDQ&ROOHJHRI,QGXVWULDO7HFKQRORJ\  <RVKLNL7$1$.$ &KXN\R0HGLFDO&RUSRUDWLRQ  <RXLFKL7$.$6+,0$ 177  
7RNLLFKLUR7$.$+$6+, 7RN\R'HQNL8QLYHUVLW\  <XNLQREX7$1,*8&+, 177  
1REXML7(7687$1, 7RN\R'HQNL8QLYHUVLW\  +LUR\XNL768-, .DQDJDZD,QVWLWXWHRI7HFKQRORJ\  +LURNR<$%86+,7$ 177  
0DVDKLUR<$1$*,+$5$ .'',5 '/DERUDWRULHV  
5\XML<$0$=$., 3DQDVRQLF&RUSRUDWLRQ  ,,((-2IILFH
2VDPX8.,*$<$
5LHNR)8.86+,0$
.\RNR+21'$

&RQWDFW,QIRUPDWLRQ
7KH,QVWLWXWHRI,PDJH(OHFWURQLFV(QJLQHHUVRI-DSDQ ,,((-㸧$UDNDZD$UDNDZDNX7RN\R-DSDQ
7HO )D[
(PDLOKHQV\X#LLHHMRUJ 
KWWSZZZLLHHMRUJ LQ-DSDQHVH 
KWWSZZZLLHHMRUJHQ LQ(QJOLVK 
KWWSZZZIDFHERRNFRP,,((- LQ-DSDQHVH 
KWWSZZZIDFHERRNFRP,,((-( LQ(QJOLVK䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷䚷,,((-7UDQVDFWLRQVRQ䚷䚷䚷䚷䚷䚷,PDJH(OHFWURQLFVDQG9LVXDO&RPSXWLQJ
䚷䚷䚷9RO1R䚷'HFHPEHU
䚷䚷&217(176

6SHFLDO,VVXHRQ-RXUQDO7UDFN3DSHUVLQ,(9&


8SRQWKH6SHFLDO,VVXHRQ-RXUQDO7UDFN3DSHUVLQ,(9&

2VDPX8&+,'$&RQWULEXWHG3DSHUV0XOWL9LHZ6WHUHR%DVHGRQ2PQLGLUHFWLRQDO(SLSRODU*HRPHWU\<XNL6$.$085$+LGHKLNR6+,6+,'2,WDUX.,7$+$5$



6\QWKHVL]LQJ+LJK5HVROXWLRQ&ORXGV8VLQJ'DQG'6LPXODWLRQV

5\RVXNH.$:$=80,<RVKLQRUL'2%$6+,

$5HJLVWUDWLRQ0HWKRGRQ3LHFHZLVH9RWLQJIRU/RZ2YHUODSSLQJ'3RLQW&ORXGV/XLV3(5$/7$0XQHWRVKL,:$.,5,.L\RVKL7$1$.$

%DFNJURXQG0XVLF6HDUFK6\VWHPWRDQ,QSXW9LGHR8VLQJ)DFWRU$QDO\VLVIRU,PSUHVVLRQ:RUGV

7RPRND]X,6+,.$:$

5HJXODU6HFWLRQ&RQWULEXWHG3DSHUV

2SWLPLVLQJ99&(QFRGLQJ8VLQJ.H\)UDPH6HOHFWLRQ

6DQDW1DJDUDMX)DQ=KDQJ6HLVKL7DNDPXUD'DYLG5%XOO$QQRXQFHPHQWV&DOOIRU3DSHUV6SHFLDO,VVXHRQWKHWK,,((-,QWHUQDWLRQDO&RQIHUHQFH
RQ,PDJH(OHFWURQLFVDQG9LVXDO&RPSXWLQJ ,(9&&DOOIRU3DSHUV6SHFLDO,VVXHRQ'DWD,QWHUZRUN6KDULQJDQG5HXVLQJ
7HFKQRORJLHV6XSSRUWLQJ'LJLWDO7UDQVIRUPDWLRQ(UD&DOOIRU3DSHUV6SHFLDO,VVXHRQ,PDJH5HODWHG7HFKQRORJ\IRU'LJLWDO)DEULFDWLRQ*XLGHIRU$XWKRUV*XLGDQFHIRU3DSHU6XEPLVVLRQ

3XEOLVKHGWZRWLPHVD\HDUE\WKH,QVWLWXWHRI,PDJH(OHFWURQLFV(QJLQHHUVRI-DSDQ 䠄,,((-䠅䠈$UDNDZD䠈$UDNDZDNX䠈7RN\R䠈-DSDQ
7HO䠖䚷)D[䠖䚷(PDLO䠖KHQV\X#LLHHMRUJ䚷KWWSZZZLLHHMRUJIIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

8SRQWKH6SHFLDO,VVXHRQ -RXUQDO7UDFN3DSHUVLQ,(9&
Editor: Osamu UCHIDA(Tokai University)
The 7th International Conference on Image Electronics and Visual Computing(IEVC2021) was held on September 8-11, 2021, as the international academic event of
Image Electronics Engineers of Japan (IIEEJ). It was based on the great success of theprevious six workshops in 2007 (Cairns, Australia), 2010 (Nice, France), 2012 (Kuching,
Malaysia), 2014 (Koh Samui, Thailand), 2017 (Da Nang, Vietnam), and in 2019 (Bali,
Indonesia). The conference aims to bring together researchers, engineers, developers, andstudents from various fields in both academia and industry for discussing the latest studies,
standards, developments, implementations, and application systems in all areas of image
electronics and visual computing. The IEVC2021 Committee initially planned to hold
IEVC2021 as a face-to-face conference in Shiretoko (Shari), Hokkaido, Japan, a World
Natural Heritage Site. However, due to the worsening of the COVID-19 infection
nationwide, the committee has decided to hold IEVC 2021 in a fully online format,
prioritizing the safety of the participants.
As well as in IEVC2019, IEVC2021 has two paper categories: general paper and latebreaking paper (LBP). Moreover, the general paper category has two tracks: Journal track(JT) and Conference track (CT). Journal track provides authors with the advantage of
having their papers published in the special issue of IIEEJ Transactions on Image
Electronics and Visual Computing scheduled for the December 2021 issue.
This special issue contains 4 papers that were accepted within the publication scheduleof December issue, and 2nd special issue on JT papers and also the special issue of “Extended
Papers Presented in IEVC2021” targeted for all presentations in IEVC2021 are planned in
the coming issues. All the papers are of high-level content and will contribute to the further
development of the field of image electronics and visual computing.
Last but not least, I would like to thank all the reviewers and editors for theircontribution to improving the quality of papers. I would also like to express my deepest
gratitude to the members of the editorial committee of IIEEJ and the staff at the IIEEJ
office for various kinds of support.

40

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）
IIEEJ

Paper

-- Special Issue on Journal Track Papers in IEVC2021 --Multi-View Stereo Based on Omnidirectional Epipolar Geometry
Yuki SAKAMURA† , Hidehiko SHISHIDO†† , Itaru KITAHARA†† (Member)† Ph.D Program in Empowerment Informatics, University of Tsukuba ,†† Center for Computational Science, University of Tsukuba<Summary> This paper introduces multi-view stereo (MVS) algorithm for an omnidirectional image(ODI). A data and smoothing term for MVS are calculated without distortion of ODIs at high latitudes
based on its epipolar geometry. Matching a normal vector with a pole with an ODI removes the distortion on
an epipolar curve. Since semi-global matching (SGM) is a one-dimensional smoothing method, a combinationof SGM can wholly remove from the distortion of the ODI. As a result of experiments, disparity and depth
can be calculated using the proposed method.
Keywords: multi-view stereo, omnidirectional image, equirectangular image, epipolar geometry, semi-global
matching

estimation method for camera model such as the essen-

1.

Introduction

tial matrix and lens distortions using the corresponding
points of omnidirectional images captured from multiple

An omnidirectional projection is prevalent in various
ﬁelds of computer vision. Mostly it is famous for AR/MR

viewpoints. This method applied Fitzgibbon’s method8)

applications which can observe all directions. Furthermore, combining an omnidirectional image (ODI) and ato the omnidirectional image captured with unknown external and internal camera parameters. Corresponding

dense depth map can realize advanced applications suchas displaying 6 degrees of freedom and reconstruction. In
the ﬁeld of stereo vision, too, various research has beendeveloped to calculate disparity using ODIs1)–5) . However, they have not entirely removed the distortion ofpoints are estimated using an epipolar geometry basedon geometric information between cameras (external parameters). A region for stereo matching is limited on theepipolar line. This epipolar constraint can shorten timecomplexity. For ODIs, the epipolar constraint is a curve
rather than a straight line9),10) .
Furthermore, epipolar geometry has another constraint
called cheirality11),12) . It reduces mismatched correspondence using an assumption which the corresponding point

ODIs. That problem is caused by a spherical projection
surface of the ODI as a ﬂat two-dimensional Markov random ﬁeld (MRF). This paper proposes a new multi-view
stereo (MVS) algorithm using omnidirectional epipolar

is not behind the camera. Cheirality can also apply to
ODIs11) . However, since the ODI has a projection surface that can capture the spherical image, epipoles always
exist. Therefore, in most MVS using ODIs, the matchingrange is limited by chirality5),10),13),14) .

geometry (OEG), which is a formulation of epipolar geometry in the equirectangular coordinates. That geometry altogether removes the distortion of the ODI by describing the neighborhood relationship of the pixels as a
one-dimensional MRF based on the spatial resolution atlatitude 0 (equator). Proposed MVS formulates a data
term and soothing term based on OEG.

2.
2.1

MVS for the conventional directional camera can simplify epipolar constraint and cheirality constraint by image rectiﬁcation. Rectiﬁcation can calculate epipolar con-

Related Worksstraints and chirality constraints without explicit calculation of epipolar geometry. For example, the omnidirectional camera can rectify two images by matching the
epipole with its poles5),14) . Also, by arranging two cameras vertically at the time of taking, rectiﬁcation can beachieved in advance, too13) .
Therefore, most of the conventional research has been

Multi-view stereo for an ODI

MVS algorithms for a conventional directional camera(such as perspective projection) use epipolar constraint6) .
In stereo vision, multiple cameras take the same environment, and the disparity is calculated using an estimatedcorresponding point of each image. Micusik7) proposes an

41

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

solved by applying the directional image to problems of
ODI. However, the distortion of the ODI has not beensolved completely. ODIs have a property of more signif-

the ground. Our proposed method is not aﬀected by the
distortion and the orientation of the ERI at all since theepipolar curve are mapped onto an equator depending on

icant distortion at higher latitude points. This problem
remains after rectiﬁcation. In He’s method15) , the ODI ismapped into a cylindrical image to eliminate the distortion. Cylindrical images are suitable to represent ODIs.

the spatial resolution of the ERI.
In Coors’s method19) , the ODI is described so thatthe sampling points are sampled evenly on the threedimensional space’s spherical surface. That method is the

However, it is not suﬃcient to describe the spatial neighborhood and correlation between pixels at high latitudes.
MVS and various computer vision ﬁelds have been
dealt with the distortion of ODIs16)–19) . The most major method is mapping an ODI to a regular polyhedron16) . Particularly in the mapping to a regular hexahedron (cube mapping), it is easy to apply the conventional image processing algorithm since the mapped planeclosest to the proposed method. The proposed method
altogether removes the distortion of the ODI by describing the neighborhood relationship of the pixels as onedimensional MRF based on the spatial resolution of thespherical surface at the equator.

becomes a quadrangle. Therefore, it is widely used in various ﬁelds17) . On the other hand, since mapping to themloses each surface’s continuity, care must be taken whenas MRF and optimizes its maximum a posteriori to estimate the disparity map. However, minimization of all
pixels’ costs is NP-hard. Therefore, it can be optimized
only when the smoothing term has submodularity.
Graph cut (GC) is often used in MVS as a highly accurate algorithm1),23),24) . Moreover, it is often used in
ODI stereo, too10),15) . However, since the smoothing

2.2

Smoothing term for multi-view stereo

MVS algorithm for directional camera treats an image

reprojecting to the ODI. Furthermore, the distortion remains near the mapped plane’s edge region since it is a
projection from a non-linear space to a linear space. In
Cheng’s method18) , six faces are connected in consideration of the connection of each face, but distortion is notterm must have submodularity in GC, the algorithm isconsidered.
Currently, the estimation of equirectangular depth

limited.
Belief propagation (BP) also is used in MVS2) . It calculates the total cost propagated from neighborhood pixels
and selects the minimum. It is also used in ODI stereo13) .
Semi-global matching (SGM) aggregates costs fromsome directions and optimizes disparity maps3),4) . As aresult, it operates faster than GC and BP. On the othermap using convolutional neural network (CNN) is active20)–22) . RectNet20) is a CNN autoencoder to estimate the depth of an equirectangular image (ERI) andis trained using 3D60 dataset20) . This research follows
Su’s method22) and applies convolution ﬁlters of diﬀerent
sizes depending on the latitude since the ERI have thehand, it is less accurate than them because the optimization is an approximate solution.
Pixels and their neighborhood relationships of the ODI
cannot be deﬁned as vertices and edges of MRF, unlike directional images, since they have the property thathigh latitude associate with signiﬁcant distortion. This

property of high distortion at high latitude. The ﬁlter
is extended horizontally at high latitude while keeping
the area constant. In contrast to MVS method, there is
no epipole on the images since the method does not use
image pairs but on-shot images. Therefore, the estimation of the disparity map converges quickly. However, the2D convolution is an approximation since the ERI has aspherical projection surface. Therefore, the edge of theﬁlter has a slight distortion. In Won’s method21) , the costproblem is caused by deﬁning a non-linear spherical projection plane as a linear two-dimensional MRF. If optimization of omnidirectional stereo can be deﬁned as a
one-dimensional optimization problem, we can remove
the distortion at high latitudes. Since SGM’s cost aggregation run independently in each direction as a onedimensional optimization, it is suitable for ODIs with
non-linear projection planes. The distortion at high latitude is removed by mapping each one-dimensional line tovolume generated using multi-view stereo is convolved for
depth estimation. The limitation of CNN method is that
it is diﬃcult to estimate the depth in untrained environments. RectNet is trained on 3D60, an indoor dataset,
which makes it diﬃcult to estimate in outdoor environments. In addition, since most of datasets consists only ofan equator. This paper proposes a new method to map
a one-dimensional line on the ODI to the equator and a
new penalty function for SGM.

images captured horizontally to the ground, they cannot
generate the depth for the ERI that are not horizontal to42

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Fig. 2 The equirectangular coordinate and the rectangular coordinate

P is the set of points with the rectangular coordinates.
Conversely, a mapping f −1 : P → M from rectangularcoordinate p to equirectangular coordinate m is given by
Eq. (3).

r= p2x + p2y + p2zpxpyθ=arctan, φ = arcsin
(3)
pz
r

Fig. 1 Epipolar geometry for the ODI. The red line isan actual limitation range with two constraints:
epipole and cheirality

3.

Omnidirectional Epipolar Geometry

3.1 Epipolar geometry for ODI
Epipolar geometry can be calculated using the cameras’ relative positions and orientation. Figure 1 shows
an overview of the epipolar geometry for the ODI. For

The epipolar constraint for the ERI is Eq. (4) using
essential matrix E and mapping f .

any pixel pb , a corresponding pixel pm is estimated. The
f (mm )T Ef (mb ) = 0

range of pm is limited in [em , qm ] by the epipolar constraint and the chieralitty constraint, where qm is an intersection point between the spherical surface of a spherical surface and straight lines parallel to the line cb p
and passes through cm , an essential matrix E is given
Eq. (1), where [R|t] is the coordinate transformation
matrix from cm to cb .
E =[t]× R

⎡

−tz
0

ty
⎥
−tx ⎦

−ty

tx

0

⎢

0

When the coordinate mb of the ERI is given, the constraint which the corresponding point mm of the otherimage must satisfy is obtained. Since the essential matrix E and the equirectangular coordinate mb are known,
they can be calculated in advance and are given by
Eq. (5).
f (mm )T N =0⎤[t]× =⎣ tz

N =Ef (mb ) = [nx , ny , nz ]T

(5)

(1)

f (mm ) is rectangular coordinate pm and deﬁnes it(pm x, pm y, pm z). Equation (6) is given from pm andEq. (5).

The ERIs is one of the projection methods for representing ODIs. In the ERIs, the horizontal axis is thelongitude, and the vertical axis is the latitude. Figure 2
shows the equirectangular coordinate m(r, θ, φ) and therectangular coordinate p(px , py , pz ).
A mapping f : M → P from equirectangular coordinatem ∈ M to rectangular coordinate p ∈ P is given by

pmx nx + pmy ny + pmz nz = 0(6)Equation (6) is the formula of the plane passing
through the origin, that is, the epipolar plane. Furthermore, N is the normal vector of the epipolar plane. Using
Eqs. (3) and (6), the epipolar curve on the ERI is given
by Eq. (7).

Eq. (2).
px =r cos φ sin θ
py =r sin φ(4)rnx cos φm sin θm + rny sin φm + rnz cos φm cos θm = (07)(2)(r, θm , φm ) are equirectangular coordinate mm mappedfrom pm . Equation (8) is obtained from Eq. (7).

pz =r cos φ cos θr is a radius of the sphere, ordinarily 1. M is theφm = arctan −set of points with the equirectangular coordinates, and43nx sin θm + nz cos θmny(8)IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Fig. 3 A normal vector of an epipolar plane and an
epipolar curve
Fig. 4 Removal for distortion of equirectangular projection. A normal vector is conformed to a pole ofequirectangular projection, especially zenith

The epipolar constraint for ERIs is calculated using
Eq. (8). Figure 3 shows the epipolar plane, its normal
vector, and epipolar curve of the ERI.

⎡

cos x
⎢
Rz (x)=⎣ sin x
0

A depth can be calculated by triangulation. Figure 1
shows a relationship between a baseline distance b, thedepth w, and the angles of a corresponding point λ. The
disparity d is deﬁned as the diﬀerence between angles of
the corresponding points in ERIs13) . The depth w is given
by Eq. (9).
b sin(λb − d)
wb =
sin d
b sin(λm + d)
wm =
sin d
d=λb − λm⎤− sin x 0⎥
cos x 0 ⎦01

A great circle is the intersection of the sphere and aplane that passes through the center point of the sphereand is one-dimensional curve on the ODI. The proposed
method can project any curves on the sphere (great circles) onto the equator using the normal vector of theplane, which includes the curves. Furthermore, the proposed method can remove the distortion of the twodimensional plane on the ODI partially. It is calculated(9)

from the spatial positional relationship between the tar3.2

Projection onto the equator

get pixel and neighborhood pixels. The normal vector of
the plane that contains the curve mapped on the equator using Eq. (10) and repeating that operation for the

Pixels and their neighborhood relations of the ODI cannot be deﬁned as vertices and edges of two-dimensional
MRF, since it has a property that high latitude associateswith signiﬁcant distortion. The proposed method calculates them directly on the sphere using a combination

number of neighborhoods can describe the neighborhood
relationship of the two-dimensional plane.
Figure 5 shows an example of the neighborhood of
3 × 3. In this case, since there are eight neighborhoods,
eight calculation is required. However, four projections
are suﬃcient since the pixel pair located at the opposite
positions are on the same curve. When the surface areaof some one-dimensional MRF. The equator’s ODI (lineon latitude 0) matches the original sphere’s spatial resolution. This means the pixels on the equator have no
distortion. Since the ODI has symmetry in every direction since it has a spherical projection surface, distortion
can be removed from the epipolar curve by mapping it toof the spherical projection surface is suﬃciently large, adiﬀerence between the spatial resolution of the ODI andthe sphere’s actual spatial resolution is almost nothing

the equator. It is achieved by rotating the to match thenormal vector N in Eq. (5) with a pole (zenith) of the
ODI. Figure 4 shows removal for the distortion.
The rotation matrix R is given by Eq. (10). θn and
φn are the equirectangular coordinates of N . Ry and Rz
are rotation matrices for the y-axis and z-axis. Therefore,
R is deﬁned as the rotation by Euler angles in the orderof yz.
ππ− φn )Ry ( − θn )22⎡⎤
cos x 0 sin x⎢⎥
Ry (x)=⎣
01
0 ⎦− sin x 0 cos x

in the near equator. Therefore, mapping the curve ina speciﬁc direction on the equator can use the mappedneighborhood pixels instead of other mapping results. In
this case, the time complexity is reduced since only one
projection is required.

4.

Multi-View Stereo

In this section, we describe an MVS algorithm using

R=Rz (

the proposed OEG. It is assumed that the geometric information of the ODI, such as positions and orientations,is known in advance. Therefore, the proposed method
estimates dense depth maps using two ODIs and their

(10)

44

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

tensor matrix derived from the gradient of the image is
given by Eq. (15).
S(p) =G(ω)ω

φ(g) becomes larger when the gradient is larger and when
the gradient direction and the aggregation direction are
equal. In the region where the texture has a large change,
Eq. (16) is based on the same assumption as conventional

each point as the spatial resolution on the sphere.
Formulation

Calculation of optimum disparity map d can be formulated as a problem of energy minimization. An energyfunction is given by Eq. (11) using data (ﬁdelity) term
Ef and smoothing term Es .

SGM: disparity may have a large change in the range
where the texture has a large change. Therefore, P is set
to a small value in those regions to deal with the sudden
disparity change.

(11)

R(p, dp , ω, dω )=P tanh(φ(g)|dp − dω |)φ(g)=(1 − g(1 − |r − g|))Data term. A local binary pattern (LBP)25) describes
each pixel as a binary pattern using the relative intensity
values of neighborhoods. Census transform26) algorithm
is based on LBP, as shown in Eq. (12).
Lp (I) =
ω∈Ω

0
2iΩ

I(p) > I(p + ω)otherwiseLp (Ib ) ⊕ Lp+dp (Im )

(12)

and neighborhood’s disparity, and the minimum is calculated. Using a quadratic curve is theoretically justiﬁed
only for correlation using the sum of squared diﬀerences(SSD).
Several post-processing ﬁxes disparity map: median ﬁltering and cross bilateral ﬁltering based on OEG. The median ﬁltering can remove small mistakes of the objects’edge. The cross bilateral ﬁltering can ﬁx the disparity(13)

map using its reference map. In this case, an RGB image is used. The cross bilateral ﬁltering is deﬁned as
Eq. (17), where σs and σc are smoothing parameters. sp

Smoothing term. SGM aggregates data term costs fromsome directions and can minimize energy function. The
proposed method can aggregate costs in spatial resolutions, the same as sphere surface using OEG. The SGM
aggregation curve is mapped on the equator, and the costs
are aggregated on that line. The smoothing term for
SGM is given by Eq. (14).
Es (d) =R(p, dp , ω, dω )(16)

4.2 Post-processing
The disparity map is estimated at the sub-pixel level
using Hirschmüller’s method3) . Then, a quadratic curve
is ﬁtted through the disparity, which has the smallest cost

Where Ω is a set of vectors for neighborhood pixels of
OEG, and iΩ is an index of Ω. The data term Ef is ahamming distance of LBP and is given by Eq. (13). ⊕ isthe exclusive disjunction (XOR) operator, and dp is thedisparity for pixel p.
Ef (d) =(15)

R is given by Eq. (16), where g is a gradient vector
derived from the structure tensor: g = λ1 x1 +λ2 x2 , λ andx are eigenvalue and eigenvector of the structure tensor.
P is the penalty. r is the SGM’s aggregation direction.

geometric information.
MVS for ODIs is calculated using the proposed OEG.
The exact resolution represents the epipolar curves forE(d) = Ef (d) + Es (d)

Ix Iy
Iy2

Ix and Iy are partial derivatives of the image I for x andy, calculated using Sobel ﬁltering. G is gaussian ﬁltering.
Sobel ﬁltering and Gaussian ﬁltering are based on OEG.

Fig. 5 An example of 3 × 3 neighborhoods of a pixel
in equirectangular projection. Each color line ismapped onto an equator, respectively

4.1

Ix2
I x Iyis a kernel for smoothing diﬀerences in coordinates, andc is a spatial kernel for smoothing diﬀerences in intensities. The cross bilateral ﬁltering can ﬁx disparity maps
without smooth edges.
Cp (I)=
s(x)=e(14)ω s(p − ω)c(I(p) − I(ω))d(ω)ω s(p − ω)c(I(p) − I(ω))
x
x
− 2σ
− 2σ
s
c

, c(x) = e

(17)

p,ω

5.

The cost function R imposes penalty depending on thediﬀerence of disparity3) . This paper proposes modifying
cost function R using structure tensor. The structure

Experiments

We conducted experiments to estimate the depth mapof a stereo pair of ERIs using the proposed method. The

45

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）(a) The virtual ERIs(a) The real ERIs(b) The ground truth depth map(b) The ground truth depth map by BLK360(c) The estimated depth map(c) The estimated depth map

Fig. 6 Virtual datasets and results: The ERIs, theground truth depth map, and the estimateddepth map

Fig. 7 Real datasets and results: The ERIs, the groundtruth depth map, and the estimated depth map

using the homogeneous transformation matrix Hs . And
then, they were transformed into the local points pd using the inverse of the homogeneous transformation matrix
Hd in scan d. Figure 7 shows examples of the generated
real ERI pair, the depth map for ground truth captured

coeﬃcient P in Eq. (16) was set to 20. The images’ sizefor experiment are 1024 × 512.
5.1

Virtual datasets

Public datasets for ODIs are primarily small-scale.

by BLK360, and the estimated depth map.

There are almost not comprehensive datasets that have

pd = Hd−1 Hs ps

ODIs, dense depth maps, and geometric information.
Therefore, we generate a virtual dataset for the experiment. Blender, a 3DCG software, generated data sets:

5.3 Outdoor datasets
We conducted an experiment using outdoor images.

ERI, depth map for ground truth, camera position andorientation.
Two virtual cameras were located in the virtual space
to capture the virtual room. The depth map for groundtruth is deﬁned as the distance from the camera to themesh. Figure 6 shows examples of the generated ERI

The same scene was captured using BLK360. Figure 8
shows examples of the generated outdoor ERI pair, thedepth map for ground truth, and the estimated depth
map.
5.4 Comparison with rectiﬁcation algorithm
We compared the results with a conventional multiview stereo algorithm using the rectiﬁcation method to
calculate epipolar constraints. The calculation of the data

pair, the depth map for ground truth, and the estimateddepth map. The depth image is visualized by linearlychanging a hue: The nearest depth is blue, and the farthest depth is red.
5.2

(18)

term and the smoothing term was the same as the proposed method (LBP and SGM). The error of each pixel
was weighted based on the spatial resolution at the equator because the spatial resolution of the ERI is smaller athigher latitudes. We used MAE (mean average absolute

Real datasets

To demonstrate the generality of the proposed method,
we conducted an experiment using real images. The same
scene was captured using BLK360 (density: 10 × 10 mm@10 m, range: 0.6-60 m, a ﬁeld of view: H360◦ × V300◦ ).
The depth image was generated by projecting all theerror) for evaluation. It is calculated from the total errorof all pixels and the number of estimated pixels. The error
was deﬁned as the sum of absolute diﬀerence (SAD) be-points onto the spherical projection surface at each camera position and orientation. Equation (18) shows thecoordinate transformation of the points. The points ps
captured in scan s was transformed into the global points

tween the estimated depth and ground truth depth maps.
Table 1 shows the results of the proposed method andthe rectiﬁcation method. We evaluated three kinds of the46

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）(a) The real ERIs

Fig. 10 The result of 3D reconstruction near the top(zenith) of the image (left: rectiﬁcation, right:
proposed method)(b) The ground truth depth map by BLK360

latitude because of the self-occlusion. Because of it, thesmoothing by the proposed SGM does not work and thedepth is not improved much. Figure 10 shows the resultof 3D reconstruction near the top of the image (zenith).
The conventional SGM aggregates from the same latitude
since it aggregates straightly on the image, and the coordinates of ERIs (x, y) are related to polar coordinate(longitude, latitude). It means that in high latitude ag-(c) The estimated depth map
Fig. 8 Outdoor datasets and results: The ERIs, theground truth depth map, and the estimateddepth map

gregation range on the sphere’s surface is smaller than
low latitude.
The error could be classiﬁed into four types (textureless, occlusion, epipole, and sudden change of disparity).
They were all complex problems to deal with using geometric methods. Textureless, occlusion, and epipole
could be improved by neighborhood pixel informationsuch as pixel values and depth values. In this experi-

Fig. 9 The result of 3D reconstruction near the image
edge (left: rectiﬁcation, right: proposed method)

ment, only simple post-processing such as median ﬁlter

region: overall pixels, high latitude pixels (greater than
66.6 degrees), and low latitude pixels (other degrees).

6.

and cross-bilateral ﬁlter was applied. However, the depth
map is improved using those methods such as peak removal, plane ﬁtting, and matching of mutual information3) .

Discussion

The depth map of the ERI could be generated using the
proposed method. The results show that the proposed

6.1

method is more accurate than the rectiﬁcation method
for all the virtual and real datasets. In particular, unlike

Application

Proposed OEG can be adapted to other MVS algorithms. In data terms for MVS, other cost functions such
as SSD, SAD, and mutual information3) can be used.
When the pixel and its neighborhood relationship can be
described as one dimension line on the sphere, we can calculate OEG since it depends on the only normal vector of
the plane, including the neighborhood’s edge. Therefore,
it can be adapted into other smoothing methods such as

the conventional rectiﬁcation method, both the data term
and the smoothing term are calculated considering thecontinuity at the ERI edges. Figure 9 shows the resultof 3D reconstruction near the image edge. The result of
the proposed method is smoother than the rectiﬁcation
method.

BP and GC if pixels and their neighborhoods’ relationship can be deﬁned as sets of one-dimensional MRF on
the sphere. However, note that GC requires submodular-

The proposed method could estimate the depth at highlatitude since OEG can calculate all pixels as the exact
spatial resolution. Table 1 shows that in a virtual dataset
the proposed method improves the depth at the high lat-ity for smoothing terms. One of the solutions is evenly
sampling on the surface of the sphere. Sampling points
can be regarded as vertices of MRF, and their line segments in the space can be regarded as edges.
A dense depth map is also used for AR/MR applica-

itude than low latitude. However, in the real dataset the
proposed method improves the depth at the low latitude
than high latitude. Figures 7 and 8 show that the images capture by BLK360 have large blank in the high47

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Table 1 Comparison with rectiﬁcation method. MAE (m) was deﬁned as the SAD between the estimateddepth and ground truth depth maps
virtual datasets
latitude
rectiﬁcation
proposed

overall
0.25510.2126high
0.13980.1162

real datasets
low
0.26540.2212overall
0.77590.5964

tions such as the 6-DOF display.

low
0.6257
0.4627overall
3.739
2.683

high
4.562
3.338

low
3.739
2.683

tern Recognition, Vol. 3, pp. 15-18 (2006).
3) H. Hirschmüller: “Stereo Processing by Semi-Global Matchingand Mutual Information”, IEEE Trans. on Pattern Analysis
and Machine Intelligence, Vol. 30, No. 2, pp. 328-341 (2008).

6.2

Limitation of omnidirectional epipolar geometry
Proposed OEG can calculate without distortion in the4) R. Spangenberg, T. Langner, S. Adfeldt, R. Rojas: “Large
Scale Semi-Global Matching on the CPU”, Proc.of IEEE Intelligent Vehicles Symposium, pp. 195-201 (2014).

high latitude of ODI. However, there are some limitations.
One is a considerable time complexity. Since the distortions of the ODI are diﬀerent from each latitude, eachpixel should be mapped. Since our proposed method
is not designed for real time processing, it takes sev-5) Z. Arican, P. Frossard: “Dense Disparity Estimation from Omnidirectional Images”, Proc. of IEEE Conference on Advanced
Video and Signal Based Surveillance, pp. 399-404 (2007).
6) Z. Zhang: “Determining the Epipolar Geometry and Its Uncertainty: A Review”, International Journal of Computer Vision,Vol. 27, No. 2, pp. 161-195 (1998).

eral hours. This is because the proposed MVS maps the
epipolar curve for each corresponding point.
Furthermore, the OEG cannot use conventional programmatic techniques for the conventional directional im-7) B. Micusik, T. Pajdla: “Estimation of Omnidirectional Camera Model from Epipolar Geometry”, Computer Vision and
Pattern Recognition, Vol. 1, pp. 485-490 (2003).
8) A. Fitzgibbon: “Simultaneous Linear Estimation of Multiple
View Geometry and Lens Distortion”, Computer Vision and
Pattern Recognition, Vol. 1, pp. 125-132 (2001).

age since the pixels’ neighborhood relationship diﬀers.
An example of the problem is dynamic programming9) R. Boutteau, X. Savatier, J.-Y. Ertaud, B. Mazari: “An Omnidirectional Stereoscopic System for Mobile Robot Navigation”,
Sensors & Transducers Journal, Special Issue on Robotic andSensors Environments, Vol. 5, pp. 3-17 (2009).

(DP) for SGM. SGM can reduce aggregation of costs using DP since a result of aggregation can use for neighborhood pixels. However, OEG cannot use the same DP
since each aggregation line is slightly diﬀerent from neighborhood pixels. Therefore, parallelization, which is dif-

10) S. Fleck, F. Busch, P. Biber, H. Andreasson, W. Strasser:”Omnidirectional 3D Modeling on a Mobile Robot Using
Graph Cuts”, Proc. of IEEE International Conference on
Robotics and Automation, pp. 1780-1766 (2005).

ferent from conventional algorithms using neighborhood
relationships, is required. Therefore, our implementation

11) R. I. Hartley: “Cheirality”, International Journal of Computer
Vision, Vol. 26, pp. 41-61 (1998).
12) T. Werner, T. Pajdla: “Cheirality in Epipolar Geometry”,
Proc. of IEEE International Conference on Computer Vision,
pp. 548-553 (2001).

of SGM is parallelized in each pixel.

7.

high
1.715
1.423

outdoor datasets

Conclusion

13) H. Kim, A. Hilton: “3D Scene Reconstruction from Multiple
Spherical Stereo Pairs”, International Journal of Computer Vision, Vol. 104, No. 1, pp. 94-116 (2013).

This paper introduced a new method to describe the
ODI and MVS algorithm. We calculated a data termand smoothing term for MVS without distortion of ODIs
at high latitudes. Matching a normal vector with a pole14) S. Pathak, A. Moro, A. Yamashita, H. Asama: “Dense 3D
Reconstruction from Two Spherical Images via Optical Flowbased Equirectangular Epipolar Rectiﬁcation”, Proc. of IEEE
International Conference on Imaging Systems and Techniques,pp. 140-145 (2016).

with an ODI removed the distortion on an epipolar curve.
Since SGM is a one-dimensional smoothing method, thecombination of SGM and OEG can wholly remove the dis-

15) L. He, C. Luo, F. Zhu, Y. Hao, J. Ou, J. Zhou: “Depth Map
Regeneration via Improved Graph Cuts Using a Novel Omnidirectional Stereo Sensor”, Proc. of Workshop on Omnidirectional Vision, pp. 1-8 (2007).

tortion. As a result of experiments, disparity and depth
could be calculated using the proposed method.

16) Y. Lee, J. Jeong, J. Yun, W. Cho, K.-J. Yoon: “SpherePHD:
Applying CNNs on a Spherical PolyHeDron Representation of
360 ° Images”, Proc. of IEEE Conference on Computer Vision
and Pattern Recognition, pp. 9181-9189 (2019).

References
1) T. Taniai, Y. Matsushita, T. Naemura: “Graph Cut BasedContinuous Stereo Matching Using Locally Shared Labels”,
Proc. of the IEEE Conference on Computer Vision and Pattern
Recognition, pp. 1613-1620 (2014).

17) R. Monroy, S. Lutz, T. Chalasani, A. Smolic: “Salnet360:
Saliency Maps for Omnidirectional Images with CNN”, Signal
Processing: Image Communication, Vol. 69. pp. 26-34 (2018).
18) H.-T. Cheng, C.-H. Chao, J.-D. Dong, H.-K. Wen, T.-L. Liu,
M. Sun: ”Cube Padding for Weakly-Supervised Saliency Prediction in 360 Videos”, Proc. of IEEE/CVF Computer Vision

2) A. Klaus, M. Sormann, K.F. Karner: “Segment-Based Stereo
Matching Using Belief Propagation and a Self-Adapting Dissimilarity Measure”, Proc. of International Conference of Pat-48

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Yuki SAKAMURA

and Pattern Recognition, pp. 1420-1429 (2018).

He received his B.E. and M.E. degree in Engineering from University of Tsukuba, Japan,
in 2017 and 2019, respectively. Since 2019, hehas been a doctoral student in the empowerment informatics in University of Tsukuba.

19) B. Coors, A. P. Condurache, A. Geiger: “SphereNet: Learning Spherical Representations for Detection and Classiﬁcationin Omnidirectional Images”, Proc. of European Conference on
Computer Vision, pp. 518-533 (2018).
20) N. Zioulis, A. Karakottas, D. Zrpalas, P. Daras: “OmniDepth:
Dense Depth Estimation for Indoors Spherical Panoramas”,
Proc. of European Conference on Computer Vision, pp. 448465 (2018).
21) C. Won, J. Ryu, J. Lim: “OmniMVS: End-to-End Learningfor Omnidirectional Stereo Matching”, Proc. of IEEE/CVF
International Conference on Computer Vision, pp. 8987-8996,(2019).

Hidehiko SHISHIDO
He received his M.E. and Ph.D. degrees in Engineering from University of Tsukuba, Japan,
in 2013 and 2016, respectively. In 2016, heserved as a researcher at the Japan Instituteof Sports Sciences. Since 2017, he has been an
Assistant Professor at University of Tsukuba.

22) Y.-C. Su, K. Crauman: “Learning Spherical Convolution for
Fast Features from 360 Imagery”, Advances in Neural Information Processing Systems, pp. 529–539 (2017).
23) M. Bleyer, M. Gelautz: “Graph-Cut-Based Stereo Matching
Using Image Segmentation with Symmetrical Treatment of Occlusions”, Signal Processing: Image Communication, Vol. 22,
No. 2, pp. 127-143 (2007).

Itaru KITAHARA (Member)He received his B.E. and M.E. degrees in Science Engineering from University of Tsukuba,
Japan in 1994 and 1996, respectively. In
1996, he joined Sharp Corporation. 20002003, he was a research associate of University of Tsukuba. He received his Ph.D. in
2003. 2003-2005, he was a researcher at ATR.
2005-2019, he was an assistant professor andassociate professor at University of Tsukuba.
Since 2019, he has been a professor at theCenter for Computational Sciences, University of Tsukuba.

24) O. Woodford, P. Torr, I. Reid, A. Fitzgibbon: “Global
Stereo Reconstruction Under Second-Order Smoothness Priors”, Proc. of IEEE Conference on Computer Vision and Pattern Recognition, pp. 1-8, (2008).
25) T. Ojala, M. Pietikäinen, T. Maenpaa: “Multiresolution GrayScale and Rotation Invariant Texture Classiﬁcation with Local
Binary Patterns”, IEEE Trans. on Pattern Analysis and Machine Intelligence, Vol. 24, No. 7, pp. 971-987 (2002).
26) R. Zabih, J. Woodﬁll: “Non-Parametric Local Transforms for
Computing Visual Correspondence”, Proc. of European Conference on Computer Vision, pp. 151-158 (1994).

(Received May 28, 2021)(Revised October 15, 2021)

49

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）
IIEEJ

Paper

-- Special Issue on Journal Track Papers in IEVC2021 --Synthesizing High-Resolution Clouds Using 2D and 3D Simulations
Ryosuke KAWAZUMI† , Yoshinori DOBASHI† (Member)† Hokkaido University<Summary> To create realistic images of outdoor scenes using computer graphics techniques, synthetic clouds with high-resolution details are required. However, the computation cost for generating highresolution clouds is very expensive and time-consuming. In this paper, we address this problem and propose
an eﬃcient method to create high-resolution images of clouds. The method combines low resolution 3D
simulation with high resolution 2D simulation. The low resolution 3D simulation is fast and is used for
determining the overall shape of the clouds. Then the high resolution 2D simulation uses the result of the3D simulation and adds small-scale details to them. We demonstrate that our method achieves four to tentimes faster computation to generate realistic high-resolution images of clouds.
Keywords: clouds, ﬂuid simulation, high-resolution method

1.

to the viewing direction and the 2D cloud simulation is
executed for each of the planes. The 2D simulation is controlled so that the shape and the motion become similar
to the input 3D clouds.

Introduction

Recently, computer graphics has been used in various
ﬁelds including entertainment and architectural designs.

In the remainder of this paper, we ﬁrst discuss the re-

For some applications such as ﬂight simulators or visualassessment of outdoor scenes, clouds play an important

lated work in section 2. Section 3 explains the governing
equations of the cloud formation processes. An overview
of our method is described in section 4. Next, section 5role to enhance the realism of synthetic images. Moreover, synthesizing high-resolution synthetic clouds is required since the resolution of the display devices is in-demonstrates the usefulness of our method with some experimental results. Section 6 discusses some limitations

creasing. For movies and video games, it is often re-in our method. Finally, section 7 concludes this paper.

quested to synthesize clouds with the desired shape andmotion. Fluid analysis has become a fundamental anduseful tool to generate realistic ﬂuid phenomena including
clouds. Miyazaki et al. developed a method for generating realistic clouds based on the numerical analysis of theatmospheric ﬂuid dynamics 1). However, this method requires an expensive computational cost particularly when
we want to generate high-resolution clouds. Generating
the desired shape of clouds is also diﬃcult by using this2.

Related Work

There are many methods for creating synthetic clouds.
They are roughly classiﬁed into two categories. One uses
procedural techniques and the other is based on the physical simulation of the cloud formation process.
Methods using procedural techniques generate the density distribution of clouds, using the idea of fractals.
Methods in this category are well described in the bookby Ebert et al2).

method.
Our purpose in this paper is to address the above problem and to generate high resolution images of clouds ef-

Methods for modeling clouds without generating
three-dimensional density distributions were also proposed 3)−6) . Although the computational cost for theprocedural technique is very low and it is possible to gen-ﬁciently by combining low-resolution 3D simulation andhigh-resolution 2D simulation. First, the user designs
the overall motion and shape of the clouds by using the
low-resolution 3D simulation eﬃciently. The resulting 3D
clouds is then used as input to the subsequent 2D simulation, which generates the high-resolution turbulent de-erate the desired cloud shapes, a trial and error process
is often required to generate realistic clouds. Creating
realistic animations of clouds is also a diﬃcult task.
Methods based on physical simulation can generate

tails. We generate a set of virtual planes perpendicular

50

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

clouds by numerical simulation1),7)–10). Although thesemethods can generate realistic clouds, their computational cost is generally expensive particularly for synthe-

we can increase the eﬃciency by simulating clouds in only
the regions where clouds exist.

3.

sizing high-resolution results.
Several methods have been proposed for enhancing ﬂow
detail as a post-process. Some of them synthesize plausible turbulent motion for a low-resolution simulation using
wavelet noise11)or curl noise12),13). In recent years, existing high-resolution ﬂow data has been used for detail enhancement. Chu et al.used a Convolutional Neural Net-

Governing Equations

Before describing our method, this section explains
the numerical simulation of cumuliform cloud formation.
First, the air is heated by the ground and this creates
ascending air currents. The temperature of the rising air
currents decreases due to adiabatic cooling, so vapor in
the air parcel causes a phase transition, condenses, andthe cloud is generated. At that time, the latent heat isliberated, which creates additional buoyancy forces andwork to learn a relation between low- and high-resolution
ﬂows14). Sato et al. introduced style transfer for turbulence, where details are transferred to low-resolution ﬂows
using patch-based texture synthesis15). These methods
could be applied to the cloud simulation. Although theseapproaches synthesize plausible details, the use of noisepromotes further growth of the clouds. These processes
are described by the following partial diﬀerential equations.
In our method, we assume that the density of air is constant and the atmosphere is an incompressible and inviscid ﬂuid media. Although these are not true for the real
atmosphere, we can simulate improtant visual eﬀects under these assumptions. The motion of the atmosphere isfunctions or other data results in ﬂows that are less realistic than those from direct, high-resolution simulations.
Rasmussen proposed a method using 2D simulators on
multiple cross sections of 3D volume for simulating largescale smoke and explosions16). Since this method uses 2D
simulation only, very high-resolution results can be gen-then expressed by the following Navier-Stokes equations.
1∂u=−(u · ∇)u − ∇p + B + f ,∂tρ∇ · u=0,

erated. The motion of ﬂuids between the cross sections
is approximately generated by interpolating the motion
computed by the 2D simulators. Therefore, the eﬀect of
the 3D motion produced by this method is limited and

(1)
(2)

where ρ is the density, B is the buoyancy force, and f is anexternal forces such as wind. We assume that the groundthe results are less realistic than that obtained by direct
high resolution 3D simulation.
Horvath et al. addressed the above problem and pro-is on the horizontal xy plane and z axis corresponds tothe vertial upward direction.Then, the buoyancy force B
is expressed by the following equation.

posed a method combining 2D and 3D simulators for synthesizing high-resolution ﬂames17). This method can improve the realism by ﬁrstly computing the overall 3D motion by the 3D simulator and then producing the high resolution details with the 2D simulator. The 3D motion of
the ﬂuid is simulated eﬃciently by using a particle-based
simulator. The result of the particle-based simulator is

B = kb

T − Tamb (h)z − gqc z,
Tamb (h)(3)

then used for the subsequent 2D simulators to create ﬁnal images. A set of virtual planes perpendicular to theviewing direction is generated and a high-resolution 2D

where T is the temperature, Tamb (h) is the ambient temperature, z = (0, 0, 1) which is a vector representing the
upward vertical direction, kb is the buoyancy coeﬃcient,
g is the coeﬃcient for the gravitational force, and qc is thedensity of water droplets or clouds Tamb (h) is a functionof height and we assume that it decreases in proportion to
the height from the ground h.The ﬁrst term on the right

simulation is used for each of the planes to compute thedetailed motion of the ﬂuid. However, this method focuses on the simulation of ﬂames and is not applicable tocloud simulation.
Our method is an extension of Horvath’s method tothe simulation of clouds. Since the cloud dynamics are

hand side of Eq. (3) indicates that the buoyancy force is
proportional to the diﬀerence between the temperaturesof the rising air parcel and the surrounding air. The second term indicates the drag force due to the falling water
droplets (i.e. clouds).
Next, the phase transition between the water vapor andmuch more complicated than that of the ﬂame, we cannot
simply apply Horvath’s method to clouds. Furthermore,

the clouds is represented by the following equations.
∂qc=−(u · ∇)qc + Cc ,
∂t

51

(4)

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）∂qv=−(u · ∇)qv + Cc + Sv ,∂t
Cc =α(qv − qs )lation time and creating high resolution images of clouds
with the desired shape. Our method combines low resolution 3D simulation with high resolution 2D simulationto achieve this goal. We use Miyazaki’s method1)for boththe 3D and 2D simulations.
First, as a preprocess, the 3D simulation is used to gen-

(5)
(6)

where qc and qv are densities of clouds and water vapor
respectively,Cc is the amount of clouds generated by thephase transition, Sv is the amount of water vapor supplied
from vapor sources such as a lake, and α is the phasetransition ratio. When Cc is less than zero, the amount ofclouds qc is reduced as shown in Eq. (4). This indicatesthe evaporation of water droplets (i.e. clouds). qs is thesaturation vapor content and the following function isused18).
qs = max(A exp(−B/(T + C)), qv + qc ),erate the overall motion and shape of the clouds. Since
we use a low resolution grid for the 3D simulation, thecomputation is fast enough for the user to repeat thesimulation with diﬀerent parameter settings until the desired result is obtained. We assume that the camera, fromwhich the clouds are rendered, is ﬁxed. Once the desired
result is generated, our method stores the distributions oftemperature, cloud density, velocity, and density of wa-(7)where A, B, and C are parameters for determining thesaturation vapor density.Finally, the temperature is expressed by the following equation.
∂T= −(u · ∇)T − Γd w + QCc + ST ,∂t

ter vapor obtained by the simulation. We also computea bounding box of the clouds for each frame and store
the largest bounding box among them to be used for thesubsequent 2D simulation.
Next, our method generates a set of virtual planes that

(8)

where Γd is the adiabatic lapse rate, w is the z component of the ﬂuid velocity, and Q is the coeﬃcient of thelatent heat. ST is the distribution of temperature on the
ground. The frist term on the right in the above equation

are perpendicular to the viewing direction as shown inindicates the advection by the ﬂow ﬁeld. The temperature of the atmosphere linearly decreases in proportion

eﬃciency, the 2D simulation is executed inside the crosssection between the virtual plane and the bounding box

to the height from the ground. The second term simulates this property; the temperature of the air parcel

computed in the preprocess for the 3D simulation. We use
the temperature, cloud density, velocity, and water vapor
density obtained by the 3D simulation for the boundary

Fig. 1(b). The high resolution 2D grid is used for eachof the planes in order to add high-frequency turbulent
details to the clouds obtained by the 3D simulation. For

decreases in proportion to the vertical movement of theparcel, w, and Γ is a proportional coeﬃcient called theadiabatic lapse rate. The third term indicates the latentcondition of the 2D simulation. This keeps the consistency between the 2D and 3D simulations. We further
use the cloud density obtained by the 3D simulation tocontrol the 2D simulation so that the resulting clouds

heat generated by the phase transition. The latent heatcorresponds to heat energy to convert from water to water vapor, which is called heat of evaporation. The heatof evaporation is deﬁned as the quantity of heat required
at a speciﬁed temperature to convert unit mass of liquid into vapor. Thus, the latent heat is proportional to

become similar to those obtained by the 3D simulation.
After the 2D simulation, we obtain a set of high resolution 2D density distributions on the virtual planes. To
render the ﬁnal image, we generate 3D density distribution from these 2D distributions by linearly interpolating
the density between the virtual planes. The ﬁnal image

the amount of clouds generated by the phase transition
as shown by the third term. For more detailed descriptions, please refer to the book by Robert A. Houze19).
We numerically solve the above partial diﬀerential equations with a ﬁnite diﬀerence method except the advection
terms (the ﬁrst terms on the right in Eqs. (1), (4), (5),is rendered by taking into account multiple scattering oflight computed with the 3D distribution. In the following, we explain the details of the 3D and 2D simulations.
4.1

and (8)). For these terms, we use the semi-Lagrangian
scheme to improve the stability of the simulation1),20).

4.

Low-resolution 3D simulation

For the low-resolution 3D simulation, We numericallysolve the governing equations described in section 3. We
use the method proposed by Miyazaki et al1). The simulation space is subdivided into a low resolution 3D grid.
At each grid point, the velocity, the pressure, the vapor

The Proposed Method

An overview of our method is illustrated by Fig. 1.
The goal of the proposed method is to reduce the calcu-

52

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）(a) low resolution 3D simulation(b) high resolution 2D simulation on
cross sections(c) final image

Fig. 1 Overview of the proposed method; (a) overall shape of the clouds designedwith low resolution 3D simulation, (b) synthesis of high resolution details
by high resolution 2D simulation on the virtual planes perpendicular to theview direction, (c) the ﬁnal image

density, the cloud (water droplets) density, and the temperature are assigned.
The initial status of the simulation is speciﬁed as fol-

3D simulation
2D simulation inside
only this regionlows. The initial temperature is set to the ambient temperature Tamb . The initial water vapor is set to decrease
exponentially from the bottom of the simulation space.
The initial temperature and water vapor are constant in
the horizontal direction. Note that this assumption is

copy density,
temperature, vapor,velocity at boundary

used for the initial conditions only and, once the simulation starts, the temperature and the water vapor are nolonger constant in the horizontal direction. The initial

Fig. 2 High resolution 2D simulation with the temperature, cloud density, vapor and velocity at theboundary, copied from the result of the 3D simulation

cloud density is zero. The initial velocity is also set tozero. A periodic boundary condition is used in the horizontal direction and a ﬁxed boundary condition is usedon the bottom and top of the simulation space.

resolution 2D grid. In order to generate similar clouds tothose obtained by the 3D simulation, the boundary con-

The user can adjust the shape and motion of the synthetic clouds by modifying the distributions of water vapor and temperature, Sv and ST , on the ground.
4.2

dition of the 2D simulation is determined by using theresult of the 3D simulation. That is, the temperature,cloud density, and water vapor density at the boundaryof the 2D grid are ﬁxed by those obtained by the 3D

High-resolution 2D simulation

The turbulent details are synthesized by the the 2D
simulation for each of the virtual planes as we mentioned
in section 4.
Since our goal is to create synthetic images of clouds,

simulation at each time step of the simulation. The turbulence is increased by using the vorticity conﬁnement
forces together with a turbulent noise function17).
We would like to mention that using the 3D simulation
data for only the boundary of the 2D simulation produces
the high-resolution turbulent details. The low-resolution

we do not have to simulate the cloud formation process
for the regions where no clouds are generated. So, we
ﬁrst identify the region where clouds are generated by

3D simulation contains only the low-frequency component of the ﬂow. Therefore, if we simply transfer the 3D
simulation data into the whole region covered by the 2D

using the result of the 3D simulation. Then, we run the2D simulation only for those regions. The details aredescribed below.
We ﬁrst compute a 3D bounding box of the clouds bya thresholding operation to the density distribution computed by the 3D simulation (see Fig. 2). We then compute a cross section of the bounding box with each of thevirtual planes. The high resolution 2D clouds are gener-

simulation, the low-frequency component is blended withthe high-frequency details obtained by the 2D simulationand the ﬁnal results are fairly blurred. We avoid thisproblem by using the 3D data for the boundary only.
Although the above method can synthesize the highresolution turbulent details, the resulting shape of theated inside only the cross section by generating a high-clouds could become diﬀerent from the 3D simulation.

53

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Fig. 3 A snapshot of the controlled 2D simulation; positive (red) and negative (green) vapor supply forthe control

Fig. 4 Low resolution 3D simulation

The reason is that we constrain the 2D simulation at only
the boundary using the data computed by 3D simulation.
A straightforward solution is to increase or decrease thesimulation is executed to generate high resolution details
so that the shape and the motion are similar to those obtained by the precomputed 3D simulation. We investigate

cloud density according to the diﬀerence between 2D and3D simulations. However, this is equivalent to the blending of cloud densities obtained by the 2D and 3D simulations, resulting in removing the high frequency details.

the eﬀectiveness of our method under diﬀerent parametersettings.
First, we verify the eﬀectivness of the method describedin section 4.2 (see Fig. 5). Figure 5(a) shows the result

Our solution is to implicitly control the 2D simulation bywithout the water vapor control. The shape of the re-adjusting the water vapor according to the diﬀerence inthe cloud density. We employ the proportional control
mechanism. At each grid point for the 2D simulation,sulting clouds is diﬀerent from that shown in Fig. 4 as
expected. Figure 5(b) directly adds the cloud density
according to the diﬀerence between the 2D and 3D sim-

we compute the diﬀerence in the cloud density between
2D and 3D simulations. Our method then artiﬁcially in-ulations. Figure 5(c) uses our method; the water vapor
is added according to Eq. (9). We set β = 2 for this sim-

creases (or decreases) the water vapor in proportion to
the diﬀerence. That is, we add Rv computed by the following equation to the water vapor.

ulation. Although the shape of the clouds in Fig. 5(b)is similar to the shape of the input clouds, the turbulent small-scale details are not suﬃciently produced. Our

Rv = β(qc3D − qc2D ),method successfully generates similar clouds with suﬃ-

(9)cient turbulent details as shown in Fig. 5(c).
In our control method, the user can adjust the degree ofthe control by tweaking parameter β. We investigate theeﬀect of this parameter; we created a set of animations
with diﬀerent values ranging from 0.1 to 2.0 sampled atwhere qc3D and qc2D are the densities of clouds computedby the 3D and 2D simulations, respectively. Since theclouds are generated by the phase transition from the water vapor to the water droplet, the cloud distribution by
the 2D simulation eventually converges to that computedby the 3D simulation.

interval of 0.1. Figure 6 shows the results of this experiment (see also the supplemental video). As β increases,
the shape of the clouds becomes similar to that of the

Figure 3 shows a snapshot of the controlled 2D simulation, where sv is positive (negative) in the red (green)
region.

5.

input clouds while the suﬃcient turbulent details are successfully generated. However, the characteristic ﬂuid-likemotion tends to disappear. The choice of β is up to theuser, but we need to provide a way that allows the user
to choose an appropriate value more intuitively, which isone of our future work. We use β = 2.0 for the rest of

Results

We implemented our method on a desktop PC withan Intel(R) Core(TM) i7-3490K CPU with 32GB RAM.
We compute a hundred-frame animation. The grid sizes

the examples shown in this paper.
Next, we compare the results with diﬀerent numbersof the virtual planes used for 2D simulation; the num-

for the 3D and 2D simulations are 100 × 100 × 100 and
200×200, respectively. The images are rendered by using
a physically-based renderer, Mitsuba21)
We ﬁrst design the shape of clouds with a low resolu-

bers are 75, 37, and 18 as shown in Figs. 7(a) through7(d). For rendering these images, the density between
the virtual planes is linearly interpolated by those of thetion 3D simulation as shown in Fig. 4. Then, the 2D

54

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）(a) No control(b) Control using cloud density(c) Control using water vapor

Fig. 5 Eﬀects of our control method using water vapor for 2D simulation(a) β = 0.1(b) β = 0.4(c) β = 0.8(d) β = 1.8

Fig. 6 Results with diﬀerent values for β(a) Proposed method:N = 75 (b) Proposed method:N = 37 (c) Proposed method:N = 18(d) Full 3D smulation

Fig. 7 Clouds simulated with diﬀerent numbers of virtual planes ((a) -(c)), (d) the
result obtained by the high resolution 3D simulation(a) Input(b) Our result
Fig. 8 An example of larger clouds
Table 1 Calculation time for all frames in seconds; N :the number of virtual planes, Full3D: the highresolution 3D simuation
N = 75 N = 37 N = 18 Full3D
3D simulation
133.5
133.5
133.5
2411.82D simulation
436.8
211.6
109.9
-

neighboring planes. We also simulated the clouds witha high-resolution 3D grid ( 200×200×200) as shown in
Fig. 7(d). Although realistic high resolution details aresynthesized in Fig. 7(d), the resulting shape is diﬀerent
from that shown in Fig. 4. Our method successfully generates clouds with similar shape to Fig. 4 and with realistic turbulent details (see Figs. 7(a) and 7(b)). However,

total

570.3

345.1

243.4

2411.8when the number of the virtual planes is small, some artifacts are observed due to the interpolation of the density
between the virtual planes as shown in Fig. 7(c). The
calculation time for all frames is summarized in Table 1.
Compared with the high resolution 3D simulation, ourmethod signiﬁcantly reduces the calculation time. Our

culation time in the table does not include the rendering
time, and the rendering time was around 1000 sec for all

method achieves 4 - 10 times faster computation. The
clouds are rendered using only density values. The cal-

The number of the virtual planes is 36. Figures 8(a) and8(b) respectively show the input clouds and the clouds

the cases.
Figure 8 shows an example of larger clouds generated
by our method. The resolutions of 3D and 2D simulations
are the same as those for Figs. 8(a) and 8(b), respectively.

55

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

ring the simulation status (water vapor, cloud density,
temperature, velocity) between the adjacent animation
frames.
We found that choosing an appropriate weight for thevapor control (β) was not an easy task. We sometimes need a trial-and-error process for choosing the goodweight value. We need to develop an automatic method
to determine the optimal weight.
Finally, let us discuss the diﬀerence of our method with
the Horvrth’s method17)which our method is based on.
First, the Horvath’s method was designed for the simulation of ﬁre and is not directly applicable to the simulation

Fig. 9 High resolution clouds

generated by our method. Our method successfully synthesizes the high resolution details.
Finally, we generate very high resolution clouds usingof clouds. We make use of the property of the cloud formation process to increase the eﬃciency and synthesize
the turbulent details. First, unlike the ﬁre simulation,
the regions where clouds forms are limited and we canidentify it by using the result of 3D simulation. So, ourour method by increasing the resolution of the 2D simulation to 400 × 400. The result is shown in Fig. 9. We
create two images of the same clouds in a daytime and at
the sunset. The realistic images of the clouds with highmethod reduces the computation time by limiting the 2D
simulation space to that region only. Next, the naive extension of the Horvath’s method to the cloud simulation

details are successfully generated. The number of the virtual planes is 75 and the computation time is 1968.46sec.

6.

is transferring cloud density from 3D to 2D simulation.
However, this approach removes the high-resolution de-Discussion

We have demonstrated that Our method successfully
creates high-resolution images eﬃciently. However, there

tails as shown in Fig. 5. We avoid this problem by using water vapor to control the 2D simulation. This approach promotes the formation/extinction of clouds in aare some limitations which we discussed below.
One of the most important limitations is that the vir-

physically-plausible way via the phase transition processtual planes need to be perpendicular to the ground. Thisis due to the fact that clouds are generated by the vertical motion of the air. This also constrains the location

and achieves a realistic synthesis of the turbulent details.
Finally, our results are considered to be less realistic
than those obtained by the direct high resolution 3Dof the viewpoint; the distance between the viewpoint andthe clouds needs to be suﬃciently long. Therefore, ourmethod is not applicable to the nearby clouds, e.g., clouds
just above the observer. Our method is suitable for synthesizing images of distant clouds.
Another limitation is that the cloud density between

simulation since our method approximately simulates thecloud dynamics with the coarse 3D simulation followed by
the high resolution 2D simulation. However, our method
is much faster than the direct high resolution 3D simulation and we can design the overall motion eﬃciently withthe coarse 3D simulation.

the virtual planes needs to be interpolated. The artifacts
can be removed by increasing the number of the places,
but the computational cost also increases. Therefore, we

7.

Conclusion and Future WorkIn this paper, we have proposed the method for generating images of clouds with high resolution details byneed to develop a method that determines the minimum
number of the planes suﬃcient to render images without
the artifacts. We could also reduce the artifacts by using
some turbulent noise functions for interpolation16).
Our current method assumes the viewpoint is ﬁxed.
So, our method cannot apply to creating the animation
with a moving camera. This problem could be addressedby moving the virtual planes according to the viewpoint.

combining 3D and 2D simulations. We demonstrated that
our method produced high quality images while signiﬁcantly reducing the computation time.
The quality of our results depends on some parametersin the method, e.g., the number of virtual planes. One of
our future works is then the development of an automatic
method to determine the parameters in the future. We
would also like to accelerate the computation by excluding such regions from the simulation that are occluded

However, we need to preserve the coherency between the
animation frames. This could be achieved by transfer-

56

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

by the clouds in front of the camera. Finally, our method
can be signiﬁcantly accelerated by the GPU since the 2D
simulation for each virtual plane can be executed inde-

ics, Vol.36, No.4, Article 14 (2017).
15) S. Sato, Y. Dobashi, T. Kim, T. Nishita: “Example-Based
Turbulence Style Transfer”, ACM Trans. on Graphics, Vol.37,
No.4, Article 84 (2018).

pendently in parallel.

16) N. Rasmussen, D. Q. Nguyen, W. Geiger, R. Fedkiw: “Smoke
Simulation for Large Scale Phenomena”, ACM Trans. on
Graphics, Vol.22, No.3, pp703–707 (2003).

Acknowledgements17) C. Horvath, W. Geiger: “Directable, High-Resolution Simulation of Fire on the GPU”, ACM Trans. on Graphics, Vol.28,
No.3, Article 41 (2009).

The authors would like to thank the reviewers fortheir constructive comments to improve our paper. This
work was supported by JSPS KAKENHI Grant Numbers18) Y. Dobashi, K. Kusumoto, T. Nishita, T. Yamamoto: “Feedback Control of Cumuliform Cloud Formation Based on Computational Fluid Dynamics”, ACM Trans. on Graphics, Vol.27,
No.3, Article 94 (2008).

JP20H05954.

References

19) R. A. Houze, Cloud Dynamics, Academic Press (1993).

1) R. Miyazaki, Y. Dobashi, T. Nishita: “Simulation of Cumuliform Clouds Based on Computational Fluid Dynamics”, Proc.
of EUROGRAPHICS 2002 ShortPresentations, pp.405—410(2002).

20) J. Stam: “Stable Fluids”, Proc. of ACM SIGGRAPH1999, Annual Conference Series, pp.121—128 (1999).

2) D. S. Ebert, F. K. Musgrave, D. Peachey, K. Perlin, S. Worley, Texturing & Modeling: A Procedural Approach, Morgan
Kaufmann (2002).

(Received May 28, 2021)(Revised October 30, 2021)

21) W. Jakob,
Mitsuba
mitsuba-renderer.org.

3) A. Trembilski, A. Broßler: “Surface-Based Eﬃcient Cloud Visualization for Animation Applications”, Proc. of the 10-th International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision 2002 (WSCG 2002),
pp.453—460 (2002).

Renderer(2010),

http://www.

Ryosuke KAWAZUMIHe is a master student of Hokkaido University. He received his BE in information science in 2021 from Hokkaido University. His
research interest is in the visual simulation ofclouds.

4) A. Bouthors, F. Neyret: “Modeling Clouds Shape”, Proc. ofEurographics 2004 (short papers), (2004).
5) F. Neyret: “Qualitative Simulation of Convective Cloud Formation and Evolution”, Proc. of Eurographics Computer Animation and Simulation Workshop 1997, pp.113—124 (1997).
6) G. Y. Gardner: “Visual Simulation of Clouds”, Computer
Graphics (Proc. of SIGGRAPH 1985), Vol.19, No.3, pp.297–
304 (1985).

Yoshinori DOBASHI (Member)
He is a Professor at Hokkaido University,
Japan. His research interests center on computer graphics, including realistic image synthesis, eﬃcient rendering, and sound modeling for virtual reality applications. He received his BE, ME and PhD in Engineering in
1992, 1994, and 1997, respectively, from Hiroshima University. He worked at Hiroshima
City University from 1997 to 2000 as a research associate.

7) Y. Dobashi, K. Kaneda, H. Yamashita, T. Okita, T. Nishita:“A Simple, Eﬃcient Method for Realistic Animation of
Clouds”, Proc. of ACM SIGGRAPH 2000, Annual Conference
Series, pp.19—29 (2000).
8) J. T. Kajiya, B. P. Von Herzen: “Ray Tracing Volume Densities”, Computer Graphics (Proc. of SIGGRAPH1984), Vol.18,
No.3, pp.165–174 (1984).
9) R. Miyazaki, S. Yoshida, T. Nishita, Y. Dobashi: “A Method
for Modeling Clouds Based on Atmospheric Fluid Dynamics”,
Proc. of the 9th Paciﬁc Conference on Computer Graphics andApplications, pp.363—372 (2001).
10) M. J. Harris, W. V. Baxter, T. Scheuermann, A. Lastra: “Simulation of Cloud Dynamics on Graphics Hardware”, Proc.
of the ACM SIGGRAPH/EUROGRAPHICS Conference on
Graphics Hardware2003, pp.92—101 (2003).
11) T. Kim, N. Thurey, D. James, M. Gross: “Wavelet Turbulence
for Fluid Simulation”, ACM Trans. on Graphics, Vol.27, No.3,
Article 3 (2008).
12) R. Narain, J. Sewall, M. Carlson, M. C. Lin: “Fast Animationof Turbulence Using Energy Transport and Procedural Synthesis”, ACM Trans. on Graphics, Vol.27, No.5, Article 166(2008).
13) H. Schechter, R. Bridson:“Evolving Sub-Grid Turbulence for Smoke Animation”, Proc. of the 2008 ACM SIGGRAPH/Eurographics Symposium on Computer Animation,
pp.1—7 (2008).
14) M. Chu, N. Thuerey: “Data-Driven Synthesis of Smoke Flows
with CNN-Based Feature Descriptors”, ACM Trans. on Graph-

57

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）
IIEEJ

Paper

-- Special Issue on Journal Track Papers in IEVC2021 --A Registration Method on Piece-wise Voting for Low Overlapping
3D Point Clouds
Luis PERALTA† (Student Member) , Munetoshi IWAKIRI†† (Member) , Kiyoshi TANAKA††† (Fellow)† Graduate School of Medicine, Science and Technology, Shinshu University ,†† Department of Computer Science, National Defense Academy of Japan ,††† Academic Assembly (Institute of Engineering), Shinshu University<Summary> To date, there is a wide range of registration methods for 3D point clouds. Nevertheless,most approaches are inadequate for pairs of 3D point clouds with low overlapping ratios due to many
incorrect estimates of putative point correspondences when numerous outliers are present. In this particularscenario, points from the non-overlapping areas represent the outliers. In this work, we introduce a piecewise voting-based method that leverages the repeatability of the transformation matrices that result fromregistering point subsets of the overlapping areas. We show the performance of our approach under the eﬀect
of its ruling parameters through extensive ablation experiments to deﬁne an easy-to-follow usage guideline
and compare it to available state-of-the-art registration methods for low overlap 3D points clouds.
Keywords: 3D point clouds, registration, low overlapping ratio, voting-based

only the necessary overlapping ratio and few point clouds

1.

Introduction

when reconstructing an object or scene due to the highcost of data acquisition - having less overlap, the fewer
point clouds acquired. Notwithstanding, the literature
has shown that CRAs are inadequate in low overlapping
ratio scenarios2),6),7), proving that these methods’ alignment quality worsens as ξ decreases. Principally due to

3D point cloud registration is the process that estimates
the transformation matrix T, composed by parameters
R ∈ SO(3) and t ∈ R3 , that best aligns point correspondences from two 3D point clouds with diﬀerent reference
frames to the same one. Producing a better detailed 3D
representation of the model or scene than the original

poor data association (i.e., the estimation of point correspondences) in the presence of many outliers that belong
to not-overlapping areas.
More sophisticated and common feature-based approaches8),9) perform better data association via apipeline employing keypoints and feature descriptor es-

point clouds attempt. Being particularly useful in applications of 3D image processing, computer vision, androbotics that require digital reconstruction.
Iterative Closest Points (ICP)1) is the most well-studied
registration algorithm2), with the most variations. Thesevariants improve speed, tolerance to noise and out-

timators. However, these also do not contemplate thelow overlap problem. Furthermore, these conventionalfeature-based methods can still yield wrong correspondences and incorrect alignments in such a scenario. Principally, due to the necessity of data-dependent parameter tuning, use of point-dispersion-dependent informa-liers, basin of convergence, and stability. Point-toplane registration (N-ICP)3), Trimmed ICP (TrICP)4),and Levenberg-Marquardt ICP (LM-ICP)5) are the moststraightforward variants since these present direct improvements for the error metric, data association, and error minimization techniques. Also, along with ICP, thesevariants are considered as Conventional Registration Algorithms (CRAs).
The overlapping ratio ξ represents the proportion ofpoints in the shared surfaces between the point clouds

tion from their keypoints/descriptor estimators, and thediﬀerence of point dispersion between truly overlapping
areas.
The interest in solving the problem of low overlap 3D
point clouds registration rose until recent years and isaddressed in diﬀerent manners. Despite the factors previously mentioned, there are feature-based methods focused on the low overlap10)–12). A similar approach to(overlapping areas). It is an essential factor in the registration process since the correspondences lay within theoverlapping areas. Moreover, it is common to intend for

58

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

making descriptors, but indirectly, is to isolate from or
exploit the information that neighboring points from theoverlapping areas possess to only focus on aligning the region of interest13),14). On the other hand, although widelyused in 2D image alignment15)–18), voting-based methods
are scarce in 3D6) but have achieved similar performance
to 2D methods. Furthermore, with the growing popularity of deep learning, the availability of more extensive
datasets19)–21), and powerful GPUs in recent years, thesehave also been applied to develop robust descriptors22),23), the registration task24)–27), and the low overlap registration problem28).
This work focuses on a voting-based registration
method that does not rely on point-dispersion informa-

which describes the transformation of points {pn } into{p̄n }. Hence, getting as close as possible to {qn } andminimizing the sum of squared point-to-point distances:tion, indirectly considers the information of neighboring
points, and does not need any training or high-end hardware to run. It works based on the assumption that
if we divide one point cloud into pieces, all the pieces,proaches.

By reordering P and Q considering C, and using an index
n, we can deﬁne the point sets {pn }{qn } and reduce theproblem to estimate the rigid body transformationn = 1, . . . , |C| := NC ,p̄n = Rpn + t

arg min3.

2

qn − p̄n  .

(2)

(3)

Related WorkAs stated in section 1, registration methods for low
overlapping 3D point clouds are relatively new. In thissection, we review details of the before-mentioned ap-

3.1 Descriptor-based
Wu et al.’s approach10) is based on the assumption thatthe overlapping areas of the source cloud lay on its bor-when registered to the other point cloud, perfectly alignto their truly corresponding areas and share a common
transformation. Because their registration process is not
aﬀected by the inﬂuence of neighboring outliers, the com-ders. It samples the border points and uses their corresponding Gaussian and mean curvatures as local descriptors to deﬁne the correspondences. However, according
to the authors, its dependence on local descriptors makes
it unstable on point clouds with uneven densities.
Sun et al.11) utilized point subsets from both point
clouds to decompose them via the K-SVD29) algorithm
to build sparse features, match the sparse representation,

mon transformation represents the one that best aligns
the overlapping areas, and subsequently, the point clouds.
Then by casting votes in the parameter space, we can ﬁnd
this transformation by leveraging its repeatability.
Section 2 presents a formal deﬁnition of pairwise registration. Section 3 introduces some of the latest developments on registration for low overlap 3D point clouds.
Section 4 introduces the proposed approach in detail, andsection 5 the experiments, which consisted of an ablation
study investigating the method’s behavior under various
parameter conﬁgurations to establish the usage guideline

and retrieve the solution T from registering matched subsets by the Coherent Point Drift30). Here the subsets act
as the registration unit to avoid potential mismatchesof section 6. Plus, a comparative study showing thismethod achieves similar qualitative results to state-ofthe-art (SoTA) methods.

Formally, we deﬁne the pairwise point cloud registration int the following manner: Having two point clouds,
NQP
the source P = {pi }N
i=1 and the target Q = {qj }j=1 ,along with a set of point correspondences C = {(i, j)},of their experiments employed handcrafted datasets witheven densities. Furthermore, the dictionary of K-SVD
is determined by minimizing a high non-convex functionthat does not guarantee to ﬁnd useful features for registration31).
Prokop et al.12) propose to use line features in both
point clouds, a trial-and-error search to match the largest
number of lines, and pruning the parameter space of thetransformation by ﬁnding each one separately throughthe objective is to ﬁnd the transformation matrix T, composed by a rotation matrix R and a translation vector t,that minimizes the sum of squared errors:

individual ﬁtness functions. The approach is capable ofregistering point clouds with overlapping ratios down to
20%. Nevertheless, they tailored the approach to speciﬁc2.

when using the whole point clouds, claiming good performance and robustness to noise. Nevertheless, mostPairwise Registration Deﬁnition2

qj − Rpi − t .

arg minscanning conditions, which limits its applications. They
assume that the point clouds coordinate system is aligned
to the sensor, allowing to estimate the rotation around(1)(i,j)∈C

59

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Fig. 1

Algorithm chart of the proposed method
computational complexity. According to the authors, it

the vertical z-axis only and translation limited to the x-yplane. Moreover, straight edges strictly exist in the point
clouds, restricting them even more to urban or man-madetook nearly 18 days to process Bunny from the Stanford
3D scan repository32).

environments. Furthermore, its susceptible parameters
drastically impact the registration accuracy, making it

3.4 Deep learning-based
The recently introduced PREDATOR28) is a neural network architecture that learns to detect the overlapping

highly unstable and hard to tune.
3.2

area and focuses on it to sample learning-based fully convolutional feature descriptors22),23). It learns where to
sample the feature points instead of a new method for

Neighboring interaction-based

Registration with Supervoxel Segmentation (RwSVS)
was also developed based on best-aligned subsets andignoring neighboring points. It uses an objective function
13)

more reliable descriptors, achieving remarkable results onmodel-type and laser scan datasets. Nonetheless, deepthat allegedly describes how good is the transformationfrom a pair of aligned subsets. The approach is capable oflearning-based methods demand many datasets to train
for a speciﬁc domain, making them not viable when there

producing alignments close to the optimal for model-type
and laser scan datasets. However, because it searches forindividual subsets that best ﬁt the objective function,it also ignores the shape information that neighboring
points can provide – prone to converge to local solutions

is insuﬃcient training data for a new type of dataset. Furthermore, nowadays deep learning models for 3D point
cloud processing are considerably resource-heavy and demand vast amounts of VRAM. PREDATOR is so complex that the authors reported to use a batch size 1 in all
experiments due to memory restrictions even with topof-the-line hardware.

in point clouds with large primitives or repetitive shapes.
HMRF-ICP14) does not ignore neighboring interactions
by using a Hidden Markov Random Field for inferring
correspondences via the Expectation-Maximization algorithm within the process of ICP. This approach reports
good alignments on laser scan point clouds but also failure cases due to insuﬃcient modeling of outliers.
3.3

4.

Proposed Method

The proposed method, detailed in the following subsections, is shown in the algorithm chart of Fig.1.
4.1 Descriptive metric
An essential part of our approach is a metric that describes good T’s as those that minimize the squared error between true point correspondences and maximize the

Voting-based

Peralta et al.6) proposed a voting-based approach thatleverages the high likelihood to converge to a good alignment of ICP when all pi ∈ P have a correspondence

overlapping ratio:in the target Q. It exploits the repeatability of obtaining good transformations every time a fully corresponded
point subsets align (an idea that we borrow in this work).
However, it has drastically low applicability because of its

ψ=

M SEp
.
ξ 1+λ

(4)

Where the overlapping ratio is ξ = (NC /NP ) · 100,
60

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

and derives from the union of the objective function of
Chetverikov et al.33) to set ξ, with the Mean Squared Error with Penalty13) to measure the mean squared error
between C points. The latter deﬁned as the points pi
with a neighbor qj that lays within the threshold radius
ε of two times the cloud resolution (CR) of P.
4.2

Table 1
Model
Armadillo
Bunny
DragonMain process

H. Buddha

The main process of the proposed approach breaks
down into the following four stages:(1) Making point subsets
This stage consists of making subsets Psv =N{Pssk }SS
k=1 , by setting a value of seed radius Rseed anddividing P via Supervoxel Segmentation (SVS)34).
(2) Aligning and evaluating
In this stage, each subset Pssk is rotated around thevertical axis at a rotation angle of ra degrees and registered with Q via a CRA (core registration) to obtain the

StairsDatasets

Point clouds
P: stand 330
Q: stand 0
P: bun270
Q: bun000
P: standRight 288
Q: standRight 0
P: standRight 288
Q: standRight 0
P: stairs 3
Q: stairs 0

ξgt (%)
32.72
33.62
46.45
20.37
50.98

Size
NP = 28, 415
NQ = 28, 220
NP = 31, 701
NQ = 40, 256
NP = 24, 573
NQ = 41, 841
NP = 72, 346
NQ = 78, 056
NP = 187, 959
NQ = 181, 077

most votes vmax , which leads to the best transformation
Tb . Since vmax derives from a set of similar T’s, an averaged transformation Tavg is formed from all these, andrepresents Tb . Notwithstanding, the event of obtaining
several bins with equal vmax may occur. If so, Tavg isestimated for each ot these bins, then evaluated via ψ,and lastly the Tavg that meets ψ → min becomes Tb .

corresponding T of their alignment. The pre-registration
rotation serves to change the initial position of Pssk andimprove the basin of convergence from the core registra-

4.3

Multi-scale process

As described in the ﬁrst stage of section 4.2, we set avalue of Rseed to control the size of the formed subsets via

tion. However, because the best initial position of each
Pssk is undeﬁned, we perform several core registrations
with various rotation angles, progressively increasing ra
at an angle step as , set by dividing 360◦ by a number ofrotations Nr .
Next, the descriptive metric of Eq.(4) evaluates the re-

SVS. Henceforth, each value of Rseed produces a diﬀerent
set Psv . Since there is not Rseed value that generates the
optimal set Psv , we perform stages 1 to 3 in a multimax
to ascale approach from a maximum seed radius Rseedminminimum Rseed . The latter are estimated based on amulti-scale usage of SVS35).

sulting T on P and Q altogether. At this point, because
it is still possible to obtain wrong T’s from the alignment of unrelated surfaces, we discard these if ψ < ψthr .
Where ψthr is a false transformations threshold.
(3) Binning and voting

5.
5.1

Experimental Evaluation

Implementation and datasets

We have implemented the proposed registration
pipeline in C++ using the Point Cloud Library36) for the

Having a putative good T on hand, it is changed to aseven-dimension vector representation T7D - formed byits unit quaternion a + bi + cj + dk and the translation
vector t = (tx , ty , tz ). Then, T7D is binned into an accumulator Accv . Where the number of bins Nb determines
the bin size in Accv and how loose or rigorously are differentiated the T7D ’s within this space. From this stage
the process uses the T7D representation, but to simplify

bounding boxes, SVS, and CRAs, along with OpenMP.
Experiments were carried out on the model-type
datasets Armadillo, Bunny, Dragon, and Happy Buddha
from the Stanford 3D Scan Repository32) and a pair oflaser scan point clouds of the scene Stairs from the ASL
Datasets37). Table 1 speciﬁes the employed point clouds
with their corresponding ground-truth overlapping ratio
ξgt and sizes.

the description of the method it is referred as T.
For the voting process, the votes vi of each bin bi are
increased by ν if the corresponding T exists in Accv . Otherwise, its bi is created and its votes initialized as vi = ν.

5.2

Evaluation measures

To assess any given alignment, we compare the corresponding Tb to the ground-truth Tgt through the ro-

Where ν is a voting step that allows compensating votes
as the bin size decreases, and deﬁnes as ν = Nb /cb .
(4) Solution estimation
The last stage consists of locating the bi that has thetation error as the geodesic distance between the rota  38) , andtions: Rerror = arccos tr Rb Rgt − 1 /2
the translation error as the L-2 norm between the translations: terror = tgt − tb .
61

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

5.3 Ablation study
5.3.1 Ablation experiments setup
First, we conduct a series of ablation experiments toblack the correspondences C. Fig. 7(b) also shows thatthe terror ≈ 1 m obtained in Stairs represents a misalignment, despite good Rerror .

observe the method’s performance under various combinations of its ruling parameters (i.e., the number of bins
Nb , the number of rotations Nr , and the core registration
algorithm) and determine the most optimal triplet conﬁguration. Each experiment consists of a scenario composedof (dataset, Nb , Nr , CRA), testing each dataset with four
numbers of bins Nb = {10, 30, 60, 90}, ﬁve numbers of rotations Nr = {4, 8, 16, 32, 64}, and the four CRAs as thecore registration algorithm (400 scenarios in total). For
the inner parameters, we set λ = 2 as used by Chetverikov
et al.33), based on observations of good alignments, we setψthr = 0.1 and cb = 10.
5.3.2 Ablation results
Figure 2 to Fig.6 show the Rerror and terror resultsof all scenarios. Each pair of charts corresponds to one

5.4
5.4.1

Comparative study
Comparative experiments setup

We divide our approach into four variants by the core
registration algorithm (denoted as Ours: ICP, Ours:
N-ICP, Ours: TrICP, and Ours: LM-ICP) and compare them to SoTA methods by the evaluation measures
Rerror , terror , and processing time.
Considering how relatively new low overlap registration
methods are, we can say that all the approaches describedin section 3. are SoTA. Nevertheless, we limit this study
to not compare to Prokop et al.12) since it requires exhaustive data-speciﬁc parameter tuning. Plus, most of
our datasets do not follow their particular assumptions.
Wu et al.10) and Sun et al.11) methods’ source codes were
not publicly released and implementing them by ourselves

CRA as the core registration algorithm, and the vertical axes correspond to the evaluation measures obtained
from each (Nb , Nr ) combination. At ﬁrst glance, it is obvious how terror does not variate along with the diﬀerent
triplets (Nb , Nr , CRA), only between datasets; staying
around 0.01 m (1 mm) in all model-type datasets andcan potentially cause misleading results. Hence, these
two approaches are not compared as well. The mentioned
voting-based method6) is also not taken into consideration
due to its lengthy processing time. Lastly, we also do notcompare to PREDATOR28) because of its high computational resources and training data demand.

above 1 m in Stairs, as its scale is extensively larger.
On the other hand, Rerror results give us more reliable
hints about the performance of the method. N-ICP as the

For these reasons, we compare our four variants to
RwSVS13) and HMRF-ICP14) only, which parameters we
set as initially described by the authors. We use the me-core registration operates best in Armadillo and Bunny,because N-ICP takes advantage of the distribution of thedian values of each core registration algorithm with thediﬀerent tested (Nb , Nr ) combinations in each dataset toexpress our variants more generally.

normals on the large round surfaces of these datasets.
This phenomenon also occurs in Happy Buddha as longas the accumulator is not sensitive to discriminating T’s(cf. Rerror in Fig.5(b) at Nb = 90). In the case of Dragonbecause its shape is less uniform than the other models,there is more variation in the results, but N-ICP and
LM-ICP perform best with Nb ≥ 30. As for Stairs, N-

5.4.2

Comparative results

Following the color pattern described before, Fig. 7
shows the visual comparison between the ground-truth
and the alignments obtainded from the evaluated methods, and Fig.8 shows the Rerror , terror , and processing
time comparatives. In the latter, each series represents a

ICP tends to converge to alignments with drift causedby the large planar surfaces of this dataset. Parameter
Nr shows to be dependent on Nb and, more importantly,
the core registration algorithm. Nonetheless, Nr ≥ 16 ismore stable with 30 ≤ Nb ≤ 60.
Although results are noisy and no conﬁguration is consistently best, conﬁgurations with the least Rerror andterror for each dataset are (Armadillo, 60, 32, N-ICP),

method; the horizontal axes correspond to the datasets,
and the vertical axes to each metric in logarithmic scale.
According to Fig.8(a), HMRF-ICP outperforms the rotation estimation in Armadillo with Rerror = 0◦ andDragon with Rerror ≈ 0.01◦ , but, as seen in Fig.8(b), it
is not better than our method variants in terms of translation. However, this improved translation does not represent substantial diﬀerences between both methods inthese datasets (cf. Armadillo and Dragon in Fig. 7 (b)(Bunny, 60, 8, LM-ICP), (Dragon, 30, 4, N-ICP), (Happy
Buddha, 60, 32, N-ICP), and (Stairs, 60, 4, ICP). Shownin the second row of alignments in Fig. 7, depicting in
red the source cloud P, in blue the target cloud Q, and inand 7(d)). HMRF-ICP also performs poorly in Bunnyand Happy Buddha in terms of rotation and Stairs for

62

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）(a) ICP(b) N-ICP

Fig. 2(a) ICP(c) TrICP(d) LM-ICP

Rerror and terror results in Bunny(b) N-ICP

Fig. 4(d) LM-ICP

Rerror and terror results in Armadillo(b) N-ICP

Fig. 3(a) ICP(c) TrICP(c) TrICP

Rerror and terror results in Dragon

63(d) LM-ICPIIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）(a) ICP(b) N-ICP

Fig. 5(a) ICP(c) TrICP(d) LM-ICP

Rerror and terror results in Happy Buddha(b) N-ICP

Fig. 6(c) TrICP(d) LM-ICP

Rerror and terror results in Stairs

translation.
On the other hand, RwSVS outperforms the Stairs

our method variants and RwSVS are signiﬁcantly slower
than HMRF-ICP. The latter converges to a solution
within 1 and 100 seconds in all datasets. RwSVS takes
from 100 to 1, 000 seconds (≈ 16 minutes) to converge inthe model-type datasets and about exhausting 100, 000
seconds (≈ 28 hours) in Stairs. Our variants take a median processing time within 1, 000 and 10, 000 seconds

dataset, estimating a signiﬁcantly better translation withterror ≈ 0.01 m (10 mm). However, it also converges tolocal minima in Bunny because of the large sphere-likesurfaces (cf. Fig. 7(c)).
Although we numerically compare the median values ofour variants, we note that these estimate nearly optimal
T’s in all model-type datasets but are deﬁcient in trans-(≈ 16 minutes to ≈ 3 hours) in all datasets.
A thorough examination of the time complexity of ourapproach allows us to comprehend the reason for the sub-lation estimation for the laser scan dataset. Furthermore,
we also note that our N-ICP variant generally performs
better in the model-type datasets as these have elaborate

optimal processing time performance. Our time complex maxmin− Rseed
· SSN log (Nr ) , or in theory, a
ity is O Rseed
log-linear time complexity O (n log n). In practice, ourtime complexity is dictated by the size of the input point
clouds, as it deﬁnes the range of values for seed resolution
max
to a minimum
Rseed to evaluate from a maximum Rseed
min
Rseed . Therefore, the larger the point clouds, the longer

shapes that help leverage the normals to estimate betterrotations. This pattern does not follow in Stairs because
the dataset comprises many non-dense and large planesurfaces that can be mistaken when divided into subsets.
Regarding the processing time, as portrayed in Fig.8(c),64IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）(a) Ground-truth(b) Ours(c) RwSVS(d) HMRF-ICP

Fig. 7Comparison of the ground-truth and alignments obtainded from each method(a) Rotation error comparative(b) Translation error comparative

Fig. 8Comparative study results

the processing time. Furthermore, the number of rotations Nr directly inﬂuences the bottleneck of the methodﬁguration, and results in section 5.4.2 supplement understanding the processing time performance. Our recommendations are the following:
1. Choose the core registration algorithm based on the
shape of the point clouds since no algorithm suits
all. Though, N-ICP and LM-ICP are generally goodin stage 2, as more rotations mean more T’s to process.

6.

(c) Processing time comparative

Usage Guideline

In this section, we specify a guideline to adopt our registration approach as the results in section 5.3.2 do notpresent a clear pattern for an exemplary parameter con-starting points.
2. We suggest setting 30 ≤ Nb ≤ 60 and Nr ≥ 16 to65

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

achieve the best possible alignment, as these tend toyield good results with any CRA as core registration.
3. Consider the application to utilize the method to deﬁne the parameter Nr as it directly aﬀects the processing time.
4. This approach is better suitable for oﬄine applica-

References
1) P.J. Besl, N.D. McKay: “A Method for Registration of 3-D
Shapes”, IEEE Trans. on Patter Analysis and Machine Intelligence, Vol. 14, No. 2, pp. 239–256 (1992).
2) L. Peng, W. Ruisheng, W. Yanxia, T. Wuyong: “Evaluationof the ICP Algorithm in 3D Point Cloud Registration”, IEEE
Access, Vol. 8, pp. 68030–68045 (2020).

tions to process model-type datasets.

7.

3) Y. Chen, G. Medioni: “Object Modeling by Registration ofMultiple Range Images”, The IEEE International Conferenceon Robotics and Automation (ICRA 1991), Vol. 3, pp. 2724–
2729 (1991).

Conclusion

We proposed a piece-wise voting registration method
for low overlapping 3D point clouds. It leverages therepeatability of transformations obtained from the align-4) D. Chetverikov, D. Svirko, D. Stephanov, P. Kresk: “The
Trimmed Iterative Closest Point Algorithm”, Proc. of the IEEE
International Conference on Pattern Recognition (ICPR 2002),
Vol. 3, pp. 545–545 (2002).

ment of subsets that arise from the overlapping areas and,via a multi-scale voting process, determines the transformation that best aligns a pair of 3D point clouds.
We conducted an ablation study that comprehended

5) A. Fitzgibbon: “Robust Registration of 2D and 3D Point Sets”,
Image and Vision Computing (IVC 2002), Vol. 21, pp. 1145–
1153 (2002).
6) L. Peralta, J. Sandoval, M. Iwakiri, K. Tanaka: “A Preliminary Study on Low Overlapping Unorganized Point Clouds
Registration Using Hough Voting”, Proc. of IIEEJ Technical
Meeting, Vol. 288, pp. 8–13 (2019).

extensive experimentation that showcased its behavior invarious scenarios deﬁned by its main parameters. Although it did not clearly deﬁne the ﬁttest parameter conﬁguration, the results enabled us to establish recommen-

7) A. Makadia, A. Patterson, K. Daniilidis, “Fully Automatic
Registration of 3D Point Clouds”, Proc. of the IEEE/CVF
Conference on Computer Vision and Pattern Recognition(CVPR 2006), Vol. 1, pp. 1297–1304 (2006).

dations on setting up the parameters. Furthermore, we
also conducted a comparative study to examine the main8) R. Rusu, N. Blodow, M. Beetz: “Fast Point Feature Histograms (FPFH) for 3D Registration”, The IEEE International
Conference on Robotics and Automation (ICRA 2009), pp.
3212–3217 (2009).

variants of our approach to available SoTA methods. Results showed that our approach achieves similar results
to SoTA in model-type datasets and even better in some9) F. Tombari, S. Salti, L. Di Stefano, “Unique Signatures of Histograms for Local Surface Description”, European Conference
on Computer Vision (ECCV 2010), (2010).

cases. Nevertheless, our approach might not be adequate

10) Y. Wu, W. Wang, K. Lu, Y. Wei, Z. Chen: “A New Method for
Registration of 3D Point Sets with Low Overlapping Ratios”,
Procedia CIRP, Vol. 27, pp. 202–206 (2015).

for applications that operate large laser scan datasets or
require real-time performance.
We see three aspects in which we can extend for fu-11) L. Sun, Y. Manabe, N. Yata, “Double Sparse Representation
for Point Cloud Registration”, ITE Trans. on Media Technology and Applications, Vol. 7, No. 3, pp. 145–158 (2019).

ture work: First, improving the discrimination against
incorrect transformations since there is still a possibility12) M. Prokop, S. A. Shaikh, K. Kim, “Low Overlapping Point
Cloud Registration Using Line Features Detection”, Remote
Sensing, Vol. 12, No. 1, Article 61 (2020).

to add these types of transformations to the accumulator
because of the reliance on the nearest neighbor function todetermine the correspondences by our descriptive metric.

13) L. Peralta, J. Sandoval, M. Iwakiri, K.Tanaka: “Pairwise Registration of Low Overlapping Unorganized 3D Point Clouds
Using Supervoxel Segmentation”, IIEEJ Trans. on Image Electronics and Visual Computing (IEVC 2019), Vol. 7, No. 2, pp.
78–87 (2019).

Second, decreasing the time complexity by eliminating
the direct inﬂuence of the number of rotations and only
dividing into subsets the potential overlapping areas into
pieces. Lastly, deﬁning an automatic parameter-tuning
descriptor robust to diﬀerentiate globally primitive surfaces highly present in laser scan datasets.

14) J. Stechschulte, N. Ahmed, C. Heckman, “Robust Low-overlap
3-D Point Cloud Registration for Outlier Rejection”, The IEEE
International Conference on Robotics and Automation (ICRA2019), pp. 7143–7149 (2019).

Acknowledgments

15) G. Stockman, “Object Recognition and Localization via Pose
Clustering”, Computer Vision, Graphics and Image Processing,Vol. 40, No. 3, pp. 361–387 (1987).

We want to thank the anonymous reviewers whose comments helped us drastically improve this manuscript’s16) D.P. Huttenlocher, S. Ullman, “Recognizing Solid Objects by
Alignment with an Image”, International Journal of Computer
Vision (IJCV 1990), Vol. 5, No. 2, pp. 195-–212 (1990).
17) H. J. Wolfson, I. Rigoutsos, “Geometric Hashing: An
Overview”, IEEE Computational Science and Engineering,Vol. 4, No. 4, pp. 10–21 (1997).

quality and Hugo Monzón and Jaime Sandoval for
their invaluable assistance, observations, and suggestions
throughout the development of this work.

18) Y. C. Hecker, R. Bolle, “On Geometric Hashing and the Generalized Hough Transform”, IEEE Trans. on Systems Man andCybernetics, Vol. 24, pp. 1328–1338 (1994).

66

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）
19) A. Geiger, P. Lenz, R. Urtasun, “Are We Ready for Autonomous Driving? The KITTI Vision Benchmark Suite”,
Proc. of the IEEE/CVF Conference on Computer Vision and
Pattern Recognition (CVPR 2012), pp. 3354–3361 (2012).

32) B. Curless, M. Levoy: “A Volumetric Method for Building
Complex Models from Range Images”, Proc. of SIG-GRAPH,
ACM, pp. 303–312 (1996).
33) D. Chetverikov, D. Stephanov, P. Kresk: “Robust Euclidean
Alignment of 3D point Sets: The Trimmed Iterative Closest
Point Algorithm”, Image and Vision Computing (IVC 2005),
Vol. 23, pp. 299–309 (2005).

20) Z. Wu, S. Song, A. Khosla, F. Yu, L. Zhang, X. Tang, J.
Xiao, “3D Shapenets: A Deep Representation for Volumetric Shapes”, Proc. of the IEEE/CVF Conference on Computer
Vision and Pattern Recognition (CVPR 2015), pp. 1912–1920(2015).

34) J. Papon, A. Abramov, M. Schoeler, F. Wörgötter: “Voxel
Cloud Connectivity Segmentation - Supervoxels for Point
Clouds”, Proc. of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2013), pp. 22–27 (2013).
35) L. Peralta, J. Sandoval, M. Iwakiri, K.Tanaka: “Automatic
Seed Resolution Selection in Supervoxel Segmentation for Registration of Low Overlapping 3D Point Clouds”, The IIEEJVisual/Media Computing Conference, (2019).

21) A. Zeng, S. Song, M. Nießner, M. Fisher, J. Xiao, T.
Funkhouser, “3DMatch: Learning Local Geometric Descriptors
from RGB-D Reconstructions”, Proc. of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR
2017), pp. 199–208 (2017).
22) C. Choy, J. Park, V. Koltun, “Fully Convolutional GeometricFeatures”, IEEE/CVF International Conference on Computer
Vision (ICCV 2019), pp. 8957–8965 (2019).

36) R. Rusu, S. Cousins: “3D is Here: Point Cloud Library(PCL)”, The IEEE International Conference on Robotics andAutomation (ICRA 2011), (2011).

23) X. Bai, Z. Luo, L. Zhou, H. Fu, L. Quan, C. Tai, “D3feat:Joint Learning of Dense Detection and Description of 3D Local
Features”, Proc. of the IEEE/CVF Conference on Computer
Vision and Pattern Recognition (CVPR 2020), pp.6358–6366(2020).

37) F. Pomerleau, M. Liu, F. Colas, R. Siegwart, “Challenging
Data Sets for Point Cloud Registration Algorithms”, The International Journal of Robotics Research (IJRR 2012), Vol. 31,
No. 14, pp. 1705–1711 (2012).
38) R. Hartley, J. Trumpf, Y. Dai, H. Li, “Rotation Averaging”,
International Journal of Computer Vision (IJCV 2013), Vol.
103, No. 3, pp. 267–305 (2013).

24) Y. Aoki, H. Goforth, R. A. Srivatsan, S. Lucey, “PointnetLK:
Robust and Eﬃcient Point Cloud Registration Using Pointnet,
Proc. of the IEEE/CVF Conference on Computer Vision and
Pattern Recognition (CVPR 2019), pp. 7156–7565 (2019).

(Received May 28, 2021)(Revised October 11, 2021)

25) Y. Wang, J. Solomon, “Deep Closest Point: Learning Representations for Point Cloud Registration”, IEEE/CVF International Conference on Computer Vision (ICCV 2019), pp. 3522–
3531 (2019).

Luis PERALTA (Student Member)

26) G. Dias Pais, S. Ramalingam, V. M. Govindu, J. C. Nascimento, R. Chellappa, and P. Miraldo, “3DRegNet: A Deep
Neural Network for 3D Point Registration”, Proc. of theIEEE/CVF Conference on Computer Vision and Pattern
Recognition (CVPR 2019), pp. 7193–7203 (2019).

He is a Dr. Eng. candidate in Telecommunication Systems at Shinshu University,
Nagano, Japan. He has obtained a B.S. degree in Mechatronics Engineering from Instituto Tecnológico de Celaya, Mexico, in 2014;
and a M.Eng. degree in Electronic Information Systems from Shinshu University in 2017.
His research interests include visual perception, 3D point cloud processing, and machinelearning applied to robotic systems.

27) C. Choy, W. Dong, V. Koltun, “Deep Global Registration”,
Proc. of the IEEE/CVF Conference on Computer Vision and
Pattern Recognition (CVPR 2020), pp. 2511–2520 (2020).
28) S. Huang, Z. Gojcic, M. Usvyatsov, A. Wieser, K. Schindler,“Predator: Registration of 3D Point Clouds With Low Overlap”, Proc. of the IEEE/CVF Conference on Computer Vision
and Pattern Recognition (CVPR 2021), pp. 4267–4276 (2021).

Munetoshi IWAKIRI (Member)He received the B.S. degree in Computer Science, and the M.Eng. degree in Mathematics
and Computer Science from the National Defense Academy of Japan in 1993 and 1998, respectively. In 1999, he joined the Departmentof Computer Science of the National Defense
Academy of Japan, as a Research Associate.
In 2002, he received the Dr. Eng. degreefrom Keio University, Japan. The institution
where he became a lecturer and associate professor in 2005 and 2015, respectively. He is amember of the Information Processing Society of Japan and pursues research in multimedia processing and information security.

29) M. Aharon, M. Elad, A. M. Bruckstein, “The K-SVD: An Algorithm for Designing of Overcomplete Dictionaries for Sparse
Representation”, IEEE Trans. on Signal Processing, Vol. 54,
No. 11, pp. 4311-–4322 (2006).
30) A. Myronenko, X. Song: “Point Set Registration: Coherent
Point Drift”, IEEE Trans. on Patter Analysis and MachineIntelligence, Vol. 32, No. 12, pp. 2262–2275 (2010).
31) R. Rubinstein, A. M. Bruckstein, M. Elad, “Dictionaries forSparse Representation Modeling”, Proc. of the IEEE, Vol. 98,
No. 6, pp. 1045–1057 (2010).

67

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Kiyoshi TANAKA (Fellow)
He received his B.S and M.S. degrees in Electrical Engineering and Operations Researchfrom the National Defense Academy, Yokosuka, Japan, in 1984 and 1989, respectively.
In 1992, he received the Dr. Eng. degreefrom Keio University, Tokyo, Japan. In 1995,
he joined the Department of Electrical andElectronic Engineering, Faculty of Engineering, Shinshu University, Nagano, Japan, andcurrently, he is a full professor in the academic
assembly (Institute of Engineering) of Shinshu University. He is the director of the Center for Global Education and Collaborationof Shinshu University, and the former VicePresident of Shinshu University in charge ofinternational aﬀairs. His research interests include image and video processing, 3D point
cloud processing, information hiding, humanvisual perception, evolutionary computation,
multi-objective optimization, smart grid, andtheir applications. Currently, he is the president of IIEEJ, a fellow of IIEEJ, a memberof IEEE, IEICE, IPSJ, JSEC, and so on.

68IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）
IIEEJ

Paper

-- Special Issue on Journal Track Papers in IEVC2021 --Background Music Search System to an Input Video Using Factor Analysis
for Impression Words
Tomokazu ISHIKAWA† (Member)† Toyo University, Prometech CG research<Summary> The use of video clip sharing sites and applications has become popular in the recent years.
However, when video contributors were looking for sound sources that could be used as background music,
it took a lot of time to actually listen to a number of songs to conﬁrm if they ﬁt the video. The authorof this study believes that video contributors would improve their production eﬃciency if there is a systemthat automatically lists up suitable background music for videos. In the proposed method, the features
for the videos are formulated using color and optical ﬂow histograms. Mel frequency cepstrum coeﬃcients,
which have often been used in the recent deep learning research, are used for the music features. A functionthat converts these features to linguistic evaluation values is obtained by multiple regression analysis. The
top ﬁve music ﬁles with the highest similarity to the input video are calculated to evaluate the proposedsystem. Whether they are compatible with the respondents’ ranking is then checked. The result shows thatthe proposed method via impression words can perform a more accurate retrieval than the learning networkthat directly connects video and music features.
Keywords: music search system, feature calculation, factor analysis

1.

Introduction

5DQNVW
3URSRVHG
V\VWHP

Video posting sites and apps, such as YouTube and TikTok, have become popular in the recent years. From the5DQNQG

standpoint of video contributors, they want to increase
the number of video views often for monetization. Adding
appropriate background music (BGM) and sound eﬀects(SE) and using interesting content and eﬃcient video

Fig. 1 Overview of the proposed system.

editing technique are considered eﬀective in increasing
the number of times the posted video is played. Therefore, video contributors search for a BGM that matcheseﬃciency. Given that, a search system is examined. Figure 1 shows a schematic diagram of the proposed system,%*0'%

5DQNUG

in which music and video are verbalized from their respective features, and the problem is to search for the opti-

the mood of the shot video or create their own. To ﬁnd
the BGM that matches their videos, movie creators mightsearch from a database or a dedicated site where muchmusic has been collected. The key to this search is themeta information in languages describing music, such as“refreshing” and “serious”, which is often determined bymum combination between music and video. The music
and video impression will be decided by recruiting subjects and conducting a questionnaire survey.

2.

the subjectivity of the person who collected the music.
In the case of a general user (e.g., a video poster or someuser without knowledge of music), it may take some time

Related Works

In this study, the BGM suitable for videos was calculated using video and music feature values. In this section, we describe related works on matching video andmusic.

to reach the desired melody, even if he/she searches froma music database or a dedicated site. This study believes that a system, which automatically lists up suitablebackground music for videos, would improve production

69

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

2.1

ing images. Liem et al. proposed a novel story-driven
approach to search the soundtrack for user-generated
videos6). They leveraged the ﬁlm’s knowledge of crossmodal relevance by searching story texts.
Hochin et al. obtained the main factors of sensitivity
by the SD (Semantic Diﬀerential) method 7)and factor
analysis8). They conﬁrmed the existence of the principal
factors of sensitivity common to images, videos, and music, and developed a prototype system for searching simi-

Video feature values

No established method has yet been recognized forvideo feature values. A video is a sequence of still images; hence, still image feature values are often used tomake video feature values by averaging them in the timeseries direction. Statistics, such as color mean value, variance, and color histogram, are often used as the feature
quantities of still images. Optical ﬂow, which expresses
the amount of movement of a corresponding object froma continuous image as a vector, is also widely used in theﬁelds of pattern recognition and computer vision.
Suzuki et al. proposed the VDUNV method that calculates the movement of human-like primates from natural videos with less constraints than optical ﬂow1). The
present study proposes a method for matching one after
another in a video by considering the vector ﬁeld calcu-

lar heterogeneous media data. They considered mapping
each media data according to their principal factors. The
proposed method also uses a part of this method.
Shimizu et al. proposed a method for generating music by learning a set of video and music features9). They
used Word2vec as the video features for color distribution,
movement distribution, and subject information. As mu-

lated by the VDUNV method as a feature. Since thismethod is speciﬁc to primate movement, the proposed
method uses general-purpose Optical Flow and color information so that it can be processed no matter what issical features, they used melodic musical features, namely
the number of notes, range, and major/minor percentages, and rhythmic musical features, namely the percentage of toms/snares, number of whole notes, and percent-in the image.

age of sixteenth notes. They were also using the evalu-

2.2

ation values for sensibility words as the user impression
evaluation, and proposing the impression learning using
a SOM (self organizing map) instead of multiple regres-Music feature values

Sato et al. proposed a system for retrieving musical works based on emotions2). This system introduces
the concept of an emotional valence measure, which hassion analysis. However, questions, especially about themethod of calculating the video feature values, are re-been proposed as a way to measure the emotional valenceof musical works. The user enters an emotional value,

maining in their work. They claimed that SOM insteadof multiple regression analysis was used due to the poor
accuracy in the absence of suﬃcient data, but we col-and the system calculates the similarity with the emotional value in the database and presents the musical workwith the highest similarity. The ﬁve factors of “exaltation,” “aﬃnity,” “strength,” “lightness,” and “solemnity”
were used as emotional values, and “tempo,” “tonality,”lect the necessary data for multiple regression analysis
through questionnaires. The video feature calculation isclariﬁed in the proposed method. Furthermore, a system that ranks and presents appropriate songs from a
database is constructed instead of generating songs.

“melody,” “rhythm,” “harmony,” and “beat” were usedas musical characteristics.
The most commonly used music features in the recent
deep learning research are log-mel spectrograms3)and mel
frequency cepstrum coeﬃcients (MFCCs)4). Hizlisoy etal. proposed an approach for music emotion recognition

3.

Proposed Method

The proposed system takes a video as the input andranks and outputs background music that seems appropriate for the input video. Video and music features areverbalized into sensibility words based on the assumptionthat similar impressions of video and music are appro-

based on the convolutional long and short memory deepneural network (CLDNN) architecture5). In their work, in
addition to standard acoustic features, they use features
obtained by feeding convolutional neural network (CNN)

priate combinations when verbalized. This similarity inverbalization is presented as the appropriate backgroundmusic for the input video. Figure 2 shows the overalllayers with log Mel ﬁlter bank energy and MFCCs. In the
proposed method, we employ MFCCs as music features
to map to impression words.
A lot of research has been published in the past on

diagram of the proposed system.
First, the video feature x from the video V was computed. This calculation method must be determined inadvance. For the BGMS to be searched, the music fea-

providing music to match the impressions of still or mov-70

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

the number of pixels in the section was counted; and thenumber of image strings was discounted to obtain the
average color distribution. The average Ci of the color
histogram of rank i is calculated as follows if the input
video is divided as a sequence of still images, and there
are N images:
N

Ci =Fig. 2 Overall view of the proposed system

ture y and principal factors of impression words zmusic
was calculated in advance and registered in the database.
The function f that converts the video feature values x
to the impression word zvideo and the function g that converts the music feature values y to the impression word
zmusic were determined by multiple regression analysis using a questionnaire survey, in which multiple people wereasked to verbalize their video and music impressions.
When the system was used, the video feature xinput was
obtained from the input video Vinput , and the impression
word zvideo was obtained by the function f . For the correspondence between the impression words of the video and

1 1
Ii (x, y, t),
W H N t=1 x,ywhere, Ii (x, y, t) is a function that returns 1 if the pixel
value at coordinate (x, y) in the tth still image matches
the rank i when the input video is considered as a sequential image. The entire histogram, where W and H are thehorizontal and vertical lengths of the image, respectively,
was normalized.
The Gunner Farnebäck method of optical ﬂow was employed to obtain the motion distribution statistics10). The
Gunner Farnebäck method is known as a denser algorithm for detecting optical ﬂow than the Lucas-Kanade
method11), and it can calculate the optical ﬂow for every pixel in the image. The Gunner Farnebäck method

the impression words of the music, the data is obtainedby a questionnaire asking which music is suitable for thesets an average window size that appropriately divides
the entire image. Accordingly, for an image resolution of
W × H, W H/w2 ﬂow vectors v = (vx , vy ) can be foundbetween images when setting in a square area w × w. A

video, and the function h is calculated in advance. By
applying this function h to zvideo , zinput which is a comparable vector in the impression word space for music is
obtained. The music feature yi for each of the multiple
BGMs in the database was calculated in advance. The
linguistic sensibility word zmusic,i was calculated by thefunction g. Subsequently, the similarity between zmusic,iand zinput was calculated. The best background music forthe input video was then ranked and presented. In otherhistogram for the ﬂow vector magnitude, with the number of bins set to 6, was created along with the average
Vi of the histogram for the magnitude of the optical ﬂow
averaged according to the time series. The average Diof the histogram, which was divided into six ﬂow vector
orientations, was also made.
The vector of the video feature values x is deﬁned asfollows by aligning the averages of the above histograms:

words, the search boiled down to solving an optimization
problem of the Equation 1 form:

x = {C1 , · · · , C12 , V1 , · · · , V6 , D1 , · · · , D6 }.
arg max sim (zinput , g(yi )) ,

(3)

(1)

Si

The music feature is one of the speech features that consider human auditory characteristics and uses MFCCs,

where sim(·, ·) is a function that calculates the similarityof two vectors.
3.1

(2)

which are commonly used for calculating acoustic similarity4). The number of dimensions of the MFCCs was

Video and music features

set to 12 along with the average of the color histogramsof the video features.

Re-formulating and using the video feature values 9)
by Shimizu et al. were considered for the video features.
The average of the color and motion distributions of thevideo was utilized according to the time series.
A color histogram was created by dividing the color
distribution into bins of four for each of the three RGB

3.2

Questionnaire survey

This study conducted a questionnaire survey to ask
people to verbalize their video and music impressions.
This questionnaire survey consisted of three steps: 1)
survey of music impressions; 2) survey of video impressions; and 3) ranking of compatibility between music andchannels. That is, the input video was used as a serial
number of still images; 0 to 255 were divided into four;71

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Table 1 Factors loading matrix for music impression word pairs
Impression word pairs in music Factor 1 Factor 2 Factor 3 Factor 4 Factor 5
Beautiful —Ugly
0.9740-0.0295
0.46070.14780.2013Natural —Unnatural
0.8381-0.0547
0.57490.52420.1781
Moist—
Thirsty
0.8840-0.24590.53480.50370.1327
Grateful —Sad
0.67930.0646
0.92630.4651
0.4885Warm—
Cold
0.7198-0.03680.90030.57550.4213Bright—Dark
0.66650.12900.93690.3676
0.5598
Tense—
Relaxed-0.47440.7426-0.4614-0.69820.2607
Simple—Complex0.2447-0.39160.77650.69780.0437
Strong—Weak-0.05120.9500-0.2632-0.4896
0.1992
Large
—
Small
0.04940.9356-0.1543-0.3449
0.1110Bold—
Delicate-0.15640.88670.1866-0.0900
0.4917
Heavy
—
Light
-0.53030.4078-0.8960-0.4718-0.3879
Dynamic —Static0.22900.71560.4558-0.0796
0.7821
Clear
—
Muddy
0.9141-0.1868
0.62490.29150.1226
Hard—
Soft-0.6138
0.5596-0.5976-0.7556-0.0204
Fresh
—
Stale
0.69990.04060.70540.28710.2357
Table 2 Factors loading matrix for video impression word pairs
Impression word pairs in video
Beautiful —Ugly
Natural —Unnatural
Moist—
Thirsty
Grateful —
Sad
Warm
—
Cold
Bright
—
Dark
Tense
—
Relaxed
Simple—Complex
Strong
—
Weak
Large
—
Small
Bold
—
Delicate
Heavy
—
Light
Dynamic —Static
Clear
—
Muddy
Sharp
—
Blunt
Fresh
—
Stale

Factor 1
0.94460.82410.89370.62430.55000.4882-0.74170.3432-0.4337-0.0725-0.2994-0.4802-0.47170.9854-0.27180.7330

Factor 2-0.3622-0.4202-0.4504-0.1561-0.29360.00380.8072-0.39490.78790.30470.73090.11030.8242-0.4030
0.9728-0.0311

Factor 3
0.50400.36880.56780.93060.81220.9952-0.44180.3906-0.2869-0.0545-0.0094-0.6489-0.0487
0.5902-0.1030
0.5800

Factor 4
0.26310.2036
0.26920.19130.21540.1223-0.3938
0.2194-0.8583-0.8391-0.6222-0.5781-0.41870.3050-0.2406-0.0259

Factor 5-0.6531-0.5722-0.6715-0.4113-0.4099-0.27170.8657-0.33110.59770.2706
0.52210.39430.5925-0.65080.3114-0.4258

video. In the music and video impression surveys, 50
videos and 50 music videos collected from existing web-the created site.
This study presents the best music for the video in asites were evaluated based on impression word pairs using
the literature 2). See Table 1 and Table 2 for the impression word pairs used.
The impression word pair evaluation in this questionnaire was on a seven-point scale based on the SD
method7). For example, in the case of “Beautiful” and“Ugly”, the respondents were asked to choose 1 for“Beautiful”, 7 for “Ugly”, and an intermediate value of 4
for neither. The length of both video and music was 30
s. The video was prepared in silence. This survey was
conducted online, and a special website was prepared for

ranked order; thus, the system evaluation required thecorrect answer that humans consider suitable. For thisreason, a stage of “ranking the suitability of music andvideo” was added to the questionnaire. If each of the 50
videos and 50 pieces of music prepared for this project iscombined, we will obtain 2500 combinations, and having
all the subjects conﬁrm all combinations will take time.
Therefore, the stage of “ranking the suitability of music and video” was dispersed in this questionnaire. One
subject was asked to watch the video edited to play 50
songs in each of the six videos and choose which combination matched comfortably from the ﬁrst to the ﬁfth
place. Figure 4 shows the appearance of the “Ranking

it. The response time for each page was measured to prevent fraud considering that an online response was being
acquired. A very short response time was not counted as

the compatibility of music and video” questionnaire.

a correct response. Figure 3 shows the appearance of72

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Fig. 4 A created survey site (ranking the suitability ofmusic and video)

Fig. 3 A created survey site (music impression survey)

3.3 Matching by verbalization
The seven-point rating values of the impression wordpairs for the video and background music obtained byimpression word zvideo of the input video, zinput whichcorresponds to the impression word of the music is calculated.

the questionnaire survey were averaged among the respondents. Five main factors were extracted by factor

zinput = Czvideo

analysis, which compressed the 16 dimensions of the studied impression word pairs to ﬁve. Let zvideo be the explanatory variable for the principal factor of video andzmusic be the explanatory variable for the principal factor of music. Assuming that the equation for converting
the video feature x to the impression word zvideo is linear, we will ﬁnd each element of the following matrix A

The similarity between zinput and zmusic , which are themain factor of impression word for music, is calculated toﬁnd the best music ﬁle for the input video.

4.

female undergraduate students of information science as
the respondents. The number of responses n = 17 included the ID numbers of the top ﬁve BGMs considered
suitable from the seven-level evaluation of 16 diﬀerent
impression word pairs for each of the 50 pieces of music, seven-level evaluation of 16 diﬀerent impression word

(4)

Each element of matrix B that converts the music feature y to the impression word zmusic was also obtainedby multiple regression analysis based on data of 1) music
impressions in questionnaire survey.
zmusic = By

Results

A questionnaire survey as regards background music
and video matching was conducted, with 17 male andby multiple regression analysis using data obtained at 2)
video impressions in questionnaire survey.
zvideo = Ax

(6)pairs for each of the 50 videos, and 50 diﬀerent BGMsfor each of the six videos. First, for the BGM, the sevenpoint ratings for each impression word pair were averaged(5)

The data obtained in questionnaire survey 3 on theharmony between video and music as perceived by peo-

to obtain 50 impression word pair ratings. These data
were subjected to a factor analysis, and the factor loading matrix was calculated. The number of factors was set
to 5. The diagonal rotation by the Promax method wasused to rotate the factors.
The maximum likelihood method was used in the factor

ple will be used to calculate the correspondence between
zvideo and zmusic . This matrix C (corresponding to h in
Fig. 2) is also precomputed by multiple regression analysis in advance. By applying this transformation to the73IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

estimation. The calculations converged. The cumulative
contribution to the ﬁfth factor was 94.81%. Similarly,
for videos, the seven-point ratings for each impression

music. As a system, this study considered the determination of the explanatory variables for songs to be addedin the future and obtained a function to calculate the ex-word pair were averaged to obtain 50 impression wordpair ratings. These data were subjected to a factor analysis. The factor loading matrix was calculated with ﬁvefactors. This calculation also converged. The cumulative

planatory variables for impression words by multiple regression analysis from music features. Table 3 presents
the results of the multiple regression analysis using the
50 songs prepared herein.

contribution to the ﬁfth factor was 89.68%.
The maximum likelihood method was used in the factor estimation. The calculations converged, and Table
1 presents the results of which. The cumulative contribution to the ﬁfth factor was 94.81%. Similarly, forvideos, the seven-point ratings for each impression word

4.2

The averages of the color histogram, optical ﬂow magnitude histogram, and ﬂow vector direction histogram
were used as the video features. OpenCV, which runs on
Python and is commonly used for image processing andanalysis, was used to calculate the features in this study.
All video ﬁles were prepared as MP4 ﬁles with 1920×1080
px resolution, 29.97 FPS, and a uniform video length of
30 s. When calculating the feature values, continuous stillpair were averaged to obtain 50 impression word pair ratings. These data were subjected to a factor analysis. The
factor loading matrix was calculated with ﬁve factors.
This calculation also converged, with the results shownin Table 2. The cumulative contribution to the ﬁfth factor was 89.68%.
4.1

Function deﬁnition using video featuresand multiple regression analysis

images were extracted from the video ﬁle. One still image
was extracted for every 10 frames; therefore, 90 still images were analyzed for each video to determine the video

Function deﬁnition using music featuresand multiple regression analysis

feature values. Figure 6 shows the extracted I-frames(frames that retain the information of all pixels without
interpolation during video compression) from one of the

The MFCCs were calculated as the background music
features herein. LibROSA, which runs on Python and iscommonly used for music analysis and signal processing,
was used for the feature computation. All music ﬁles arevideos used in the experiment arranged in order from the
top left. Figure 7 depicts an example of visualizing the

30 s long and converted to .wav ﬁles with a sampling rateof 44.1 kHz. The music features were calculated using

average of the color histogram for this video.
In this study, the video is the input; thus, it is necessary

sound wave data as the input. The number of dimensionsof the MFCCs was set to 12. Figure 5 shows an example

to determine the function f to calculate the objective
variable zvideo of ﬁve-dimensional impression words fromthe video features using the previously sampled video. Inof the MFCC visualization.
The background music was kept in the database; hence,
the main factor zmusic of the impression words of ﬁve dimensions can be calculated in advance for each piece ofthe questionnaire, the impression rating values for the 50
videos were calculated. Five of which were used for thesystem evaluation. The remaining 45 videos were used forthe multiple regression analysis. The multiple regression
analysis showed that six of the 24 video features could

Fig. 5 An example of MFCC visualization

Fig. 6 An example of videos used in the experiment74

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Table 3 A matrix for calculating the explanatory variables of impression words from music features (Themean value of R, the multiple correlation coeﬃcient, was 0.579)

hhh

hhh

Music Featureshhh
hhh
Impression words
hhh

Factor 1
Factor 2
Factor 3
Factor 4
Factor 5hhhhhh
h Music Features

1

2

3

4

5

-0.0003-0.0029-0.0107-0.0025-0.0011-0.0053-0.0171
0.0435-0.0200-0.00410.00820.0039-0.0144-0.01540.01650.0045-0.0184
0.0120-0.0073
0.0076-0.0074
0.0120-0.02420.00270.0069

7

8

9

10

0.04020.0103-0.0600
0.0328-0.0266

0.02320.01230.0606-0.0003-0.0265-0.0526-0.0352
0.0072-0.01130.06070.00300.00180.0947-0.0732-0.1310

12

constant-0.0371-0.0270-0.0656
0.14880.1252

3.366110.0414-5.65311.0184-5.8727hhh
6hhh
Impression wordshh
Factor 1-0.0296
Factor 2-0.0027
Factor 3-0.0095
Factor 4
0.0185
Factor 5-0.0070hhh
hhh
Music Features
hhh
11
hhh
Impression words
hhh
Factor 1
0.0341
Factor 2
0.0258
Factor 3-0.0134
Factor 4-0.0315
Factor 5-0.0313be represented by linear combinations of other factors.
Eighteen explanatory variables were used to obtain an
equation for calculating the target variable factor, whichis the impression word.
Having determined the transformation method fromeach feature, the proposed system was then evacuated.
In the system, the video features x for the input video
M were ﬁrst calculated. Next, the function f and h wereused to convert the video features to the principal factors of the impression words for music. Subsequently, thesimilarity between the main factor and the pre-calculated
music impression words in the database was calculated,
and the top ﬁve was determined. The similarity was
based on the cosine similarity as follows:

Fig. 7 An example of the average of the color histogram

be represented by linear combinations of other factors.
Eighteen explanatory variables were used to obtain ansim(zinput , zmusic ) =equation for calculating the target variable factor, whichis the impression word. Table 4 lists the results.
4.3

zinput · zmusic|zinput | |zmusic |(7)

Table 5 shows the relevant results. The numbers in
the table indicate the music ﬁle IDs. For example, forthe video 19, the system judged the music ﬁle with ID 4

System evaluation

In this study, the video is the input; thus, it is necessary

to be the best match, followed by that with ID 14, thenthat with ID 8 in the third rank. The respondents tothe survey were also asked to rank the ﬁles, for example,
Respondent 1 (Res. 1) judged that the music ﬁles of 26,
41, 23, 12 and 9 are suitable for the input video 19 in

to determine the function f to calculate the objective
variable zvideo of ﬁve-dimensional impression words fromthe video features using the previously sampled video. In
the questionnaire, the impression rating values for the 50
videos were calculated. Five of which were used for thesystem evaluation. The remaining 45 videos were used forthis order. The background color in the table indicatesthat the respondent selected the same music ﬁle as thesystem.

the multiple regression analysis. The multiple regression
analysis showed that six of the 24 video features could75

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Table 4 A matrix for calculating the explanatory variables of impression words from video features (The meanvalue of R, the multiple correlation coeﬃcient, was 0.657)

hhh

hhh

Video veatureshhh
hhh
Impression words
hhh
Factor 1
Factor 2
Factor 3
Factor 4
Factor 5

Factor 1
Factor 2
Factor 3
Factor 4
Factor 5hhhhhhh Video veatures

1

2

3

4

5

-6.74968.9573-1.3093-0.4481-16.374668.8506-11.7455
1.71680.587621.4717-2.3585
3.1300-0.4575-0.1566-5.7218
7-5.2684
6.9916-1.0220-0.3498-12.7812-5.81637.7188-1.1283-0.3862-14.11058-5.37197.1290-1.0420-0.3567-13.0324

11.7737-15.6247
2.28390.781728.56329-2.4919
3.3070-0.4834-0.1654-6.0455

8.3493-11.08021.6196
0.554320.255510-43.1290
57.2360-8.3662-2.8634-104.6319

12

13

14

15

-75.3922
100.0521-14.6247-5.0055-182.9031
17-0.9933
1.3181-0.1927-0.0659-2.4097-44.0803
58.4985-8.5507-2.9266-106.9398
18-1.9495
2.5872-0.3782-0.1294-4.7296-91.4889
121.4138-17.7471-6.0741-221.9539
constant
40.6915-54.0013
7.89342.701698.7185-3.8072
5.0525-0.7385-0.2528-9.2364hhh
11hhh
Impression wordshh
Factor 1-31.4877
Factor 2
41.7869
Factor 3-6.1080
Factor 4-2.0905
Factor 5-76.3897
16
Factor 1
4.1879
Factor 2-5.5577
Factor 3
0.8124
Factor 4
0.2780
Factor 5
10.1599

Table 5 Ranking of background music that is considered best for the video (The numbers in thetable are the IDs of the music ﬁles)XXXtargetTable 6 Comparison of evaluation index results bysearch method
MethodFrom video to music directly
Via impression words (Proposed)

XXX

video 19

video 21

video 22selector

XXX Sys. Res. 1 Res. 2 Res. 3 Res. 4

1st
2nd
3rd
4th
5th
1st
2nd
3rd
4th
5th
1st
2nd
3rd
4th
5th

4
14
8
24
12
5
9
48
41
22
5
9
2
22
15

26
41
23
12
9
24
1
3
22
23
48
15
19
30
5

13
49
38
34
24
22
17
39
24
34
30
18
15
35
32

22
24
17
20
34
33
9
25
1
42
22
30
18
26
15

41
7
20
3
14
7
17
34
36
37
43
19
40
15
30

MAP
0.483
0.653

NDCG
0.464
0.613

background music”. MAP (Mean Average Precision)12)
was calculated to be 0.653, assuming that the target music could be retrieved if music IDs which the proposedsystem outputs are included in someone’s answer. The
NDCG (Normalized Discounted Cumulative Gain) was
calculated to be 0.613. The number of times the respondent selected the BGM was used as the degree of conformity in this calculation.
The results of calculating MAP and NDGC in comparison with the retrieval system that directly connects video
and music features are summarized in the Table 6. The
retrieval system for this comparison was based on the
results of the questionnaire obtained this time, and D
was determined by multiple regression analysis to calculate the similarities sim(Dx, y). The multiple correlation

In addition, this study presented songs that were generally in the top ﬁve of those that the respondents felt
were suitable. As a qualitative evaluation of the proposed system, we asked the questionnaire respondents to
conﬁrm the ranking of the best background music output
by the system and ask their impressions, and we got thecoeﬃcients of the video and music features obtained bymultiple regression analysis are about 0.300 on average,
which is probably one reason for the low accuracy.

comment as “it is acknowledged that some of the system
output contain a good matching between the video and76

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

5.

4) B. Logan: “Mel Frequency Cepstral Coeﬃcients for Music
Modeling”, Proc. of International Symposium on Music Information Retrieval, Vol.270, pp.1–11 (2000).

Conclusion and Future Work

This study proposed a search system that uses video
as an input and presents the most suitable music ﬁle forit. Through a questionnaire survey, the impressions of5) H. Serhat, Y. Serdar, T. Zekeriya: “Music Emotion Recognition Using Convolutional Long Short Term Memory Deep Neural Networks”, Journal of Engineering Science and Technology(2020).

images and music were quantiﬁed by evaluating the language used to describe each. The study respondents wereasked to select the top ﬁve most suitable background music for the video in the questionnaire.
The video features were formulated using color and op-6) C. C. S. Liem, A. Bazzica, A. Hanjalic: “Musesync: Standingon the Shoulders of Hollywood”, Proc. of the 20th ACM International Conference on Multimedia, MM ’12, pp.1383–1384(2012).

tical ﬂow histograms, while the MFCCs were used for themusic features. A function that converted these features
to linguistic evaluation values was obtained by multiple8) T. Hochin, K. Yamada, T. Tsuji: “Multimedia Data Access Based on the Sensitivity Factors”, Proc. of 2000 International Database Engineering and Applications Symposium(Cat. No.PR00789), pp.319–326 (2000).

regression analysis. The proposed method retrieves thebest BGM via impression words, which are hidden factors in the input video. It was shown that the proposed
method improves the retrieval accuracy more than directly associating video feature values with music feature9) Y. Shimizu, S. Kanno, T. Itoh, S. Sagayama, M. Takatsuka:“Automatic Music Component Selection Based on Impression of Input Movies”, Proc. of 2018 Nicograph International(NicoInt), pp.22–25 (2018).

7) J. G. Snider, C. E. Osgood, Semantic Diﬀerential Technique;
A Sourcebook. Aldine Publishing Company (1969).

10) G. Farnebäck: “Two-frame Motion Estimation Based on Polynomial Expansion”, Proc. of the 13th Scandinavian Conference
on Image Analysis, SCIA 03, pp.363–370 (2003).

values. The obtained results will be useful for determin-

11) B. D. Lucas, T. Kanade: “An Iterative Image Registration
Technique with an Application to Stereo Vision”, Proc. of the
7th International Joint Conference on Artiﬁcial Intelligence,
pp.674–679 (1981).

ing what kind of impression words are appropriate to add
to the current BGM database when adding metadata toit.

12) C. Buckley, E. M. Voorhees: “Evaluating Evaluation Measure
Stability”, Journal of SIGIR Forum, Vol.51, No.2, pp.235-242(2017).

The future work will focus on improving the system accuracy by collecting data from a large number of people(Received May 28, 2021)(Revised September 30, 2021)

by opening the site to the public and reconsidering the
calculation method of the feature values considering thatnot many subjects can be collected by using only ques-

Tomokazu ISHIKAWA (Member)tionnaires. If the number of data can be increased, weplan to compare the proposed method with the method

He received the B.S. degrees in Information
Science from Tokyo University of Science in
2003, and his M.S. and Ph.D in Complexity
Science and Engineering from The Universityof Tokyo in 2005, 2012, respectively. He was
hired as a system engineer by IBM Japan,
Ltd. until 2007. He is currently an Associate Professor in Information Networking forInnovation and Design (INIAD) at Toyo University. His research interests center on computer graphics, including visual simulation.

using CNN. Further more, since the videos and music
used in the questionnaire were selected by the author,
it is possible that more accurate data could have beenobtained if arbitrariness had been eliminated in the selection process. In addition to examining the method offeature value calculation, sampling that can adequately
represent the feature space is also necessary.

References
1) W. Suzuki, N. Ichinohe, T. Tani, T. Hayami, N. Miyakawa,
S. Watanabe, H. Takeichi: “Novel Method of Extracting Motion from Natural Movies”, Journal of Neuroscience Methods,
pp.291:51–60 (2017).
2) A. Sato, J. Ogawa, H. Kitakami: “An Impression-Based Retrieval System of Music Collection” Proc. of Fourth International Conference on Knowledge-Based Intelligent Engineering
Systems and Allied Technologies (Cat. No. 00TH8516), Vol.2,
pp.856–859 (2000).
3) H. Purwins, B. Li, T. Virtanen, J. Schlüter, Shuo Yiin C., T.
Sainath: “Deep Learning for Audio Signal Processing”, IEEE
Journal of Selected Topics in Signal Processing, Vol.13(2),pp.206–219 (2019).

77

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Optimising VVC Encoding Using Key Frame Selection
Sanat Nagaraju†, Fan Zhang†, Seishi Takamura (Member)††, David R. Bull††Visual Information Laboratory, University of Bristol, United Kingdom††NTT Computer and Data Science Laboratories, NTT Corporation, Japan<Summary> This paper presents a simple approach to encode video sequences by selecting an optimal Intra encodedframe (I-frame). I-frames in Versatile Video Coding (VVC) random access mode influence the bitrate required for theremaining frames in the group of pictures (GOP). Hence, we carry out a full search for the optimal frame in the first GOP and
select an I-frame with appropriate spatiotemporal information. This full search is carried out during the pre-processing stage
and hence, do not introduce large complexities into the VVC codec. We prove our technique with over 100 sequences, using
BVI-DVC and JVET sequences, to have wide range of content. Our technique is best suited for offline coding and suitablefor pre-recorded streaming videos. Applying this technique to encode JVET sequences using the VVC random access mode,
results indicate average bit rate savings (based on the Bjøntegaard Delta metric) of approximately 2.1% on the Y-channel.
Keywords: video coding, reference frames, intra frame selection, coding gains, versatile video coding

spatiotemporal redundancies by selecting the optimum
reference frame. In Section 4, experiments are reported thatcharacterise our approach using the BVI-DVC database6)
and JVET sequences7). The results demonstrate significant
bitrate savings on the BVI-DVC database while showing
potential savings of over 2% on JVET sequences under
different conditions. Finally, in Section 5, we provide ourconclusions and next steps.

1. Introduction
The COVID pandemic has shown an importance ofvideo communications – from initial medical screening toentertainment video streaming, social gaming, multi-user
online meetings and many more applications. A study
conducted by Cisco in 2018 predicts that the video traffic
accounts for 80% of the internet traffic1). This imposessignificant demands on the internet bandwidth and toaddress this, video coding standards must continuously
evolve. Versatile Video Coding2) standard also known as
ITU-T H.266 | ISO/IEC 23090-3 MPEG-I Part 3 VVC, isthe latest of these standards, providing up to 50% bitrate
reduction over its predecessor, H.265/HEVC3).
VVC builds on the HEVC standard by employing moresophisticated coding tools, one of which is an extension of
the intra coding tool to support 65 angular prediction modes
in addition to planar and DC modes. These finer intra
modes combined with a quadtree split of the picture enable
an intra coded frame to be coded efficiently4). Because these
intra coded frames (I-frames) influence the coding
efficiency of subsequent frames (P-/B-frames) particularly
in Random Access (RA) mode5), the way I-frames are
selected and encoded has a significant influence on overall
performance. The traditional method of assigning I-frames
as Long-Term Reference (LTR) frames further enables
exploitation of spatiotemporal redundancy.
This paper presents an approach that optimises theselection of an I-frame from a candidate list as a basis for
improving the overall coding efficiency of VVC. Our
approach is suitable for offline encoding of video content
of any resolution which can then be transmitted in RA mode
without any additional delay at a range of bit-rates. Our
primary objective is to verify the coding gains achieved
with a selection of a suitable I-frame.
Our paper is organised as follows: first, we review theprevious techniques related to selecting optimal reference
frames/I-frames in Section 2. Next, we introduce ourtechnique in Section 3 and describe how it exploits

2. Related Work
De-Luxán-Hernández et al. propose an Intrasubpartition (ISP) coding mode for VVC8) which divides
the intra predicted blocks into 2 or 4 sub-partitions. Thismethod, adopted in VVC at the 13th JVET meeting,provided a Bjøntegaard Delta (BD-RATE) 9), 10)
improvement of 0.27% in RA mode and 0.57% in All Intra(AI) mode. The technique implements intra coding without
consideration of how the coding block impacts subsequent
frames in the Group of Pictures (GOP) / intra period 11).
Much of the prior research into reference frame
selection introduces additional frames into the sequence.
Erabadda et al.12) propose the introduction of a virtualreference frame, formed using a statistical-based machinelearning approach. They then apply the virtual reference
frame to static-camera applications with repeated shots.
However, this technique is applied to HEVC and appends
additional frames to the sequence. Takamura et al.13)
proposed a coding method which also adds additionalframes to the sequence. This method, when applied toHM16.414) provides a 57% BD-RATE gain and is limited
to static-camera sequences.
In conventional approaches, the first frame in a GOP isused as I-frame with no consideration of spatiotemporal
redundancies. Our approach improves upon this bysearching for a superior candidate I-frame, thus maximising
the spatiotemporal information for the remainder of theintra period. We further improve over the De-LuxánHernández et al. Intra subpartition (ISP) method8) byenhancing the RA mode in our approach (RA mode

78

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Intra Period

GOP 0

GOP 1

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

Coding Order

0

5

4

6

3

8

7

9

2

12

11

13

10

15

14

16

1

21

20

22

19

24

23

25

18

28

27

29

26

31

30

32

17

Slice Type

I

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

I

POC

Fig. 1 Hierarchical B Coding structure in VVC RA mode. Illustration with GOP size = 16 and Intra period = 32

GOP 0
Original Sequence

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

Candidate frame in key frame position

6

0

1

2

3

4

5

7

8

9

10

11

12

13

14

15

16

Fig. 2 Proposed Method – move the key reference frame (6th frame as the example) to I-frame position

provides significant savings when compared to All Intramode). Finally, in our approach, we remove the need forany additional frames when compared to the method
required by Erabadda et al12) and Takamura et al13), thus
saving significant bits.

coding order of the remaining GOPs is unchanged. Basedon this rule we can reduce the number of candidates
required for full search to 16 (i.e., the GOP size).
Furthermore, for each key reference frame, the VVC codec
takes care of unit-wise encoding based on its rate distortion
optimisation process; hence, it is expected that the I-frame
along with subsequent frames benefit from this enhanced
spatiotemporal information. This is the default behaviour of
the VVC codec, and we have not altered this.
For each key reference fame candidate, as illustrated inFig. 2, we move it to the first I-frame position, then encode
the sequence, and evaluate the coding efficiency in BDRATE compared to the anchor (the unchanged sequence).
We employ the BD-RATE measure as it encompasses bothbit-rate and PSNR change. BD-RATE measurements arecommonly used by JVET and the research community forcomparing video coding performance7) as it provides a good
method of evaluating the overall rate quality difference
between two coding approaches.
The best key reference frame is chosen from the 16
candidates. The standard compliant encoded bitstream isthen transmitted to the decoder. An additional single byte isalso transmitted to signal the key frame position at the
decoder. The decoder uses the bitstream and the key frame
position to restore the original frame order. During
encoding, the cost of this approach relates to a full searchin the first GOP.
Our primary target application is offline compression
using the RA mode, as it is typically used for DVD videos
and streaming services such as YouTube and Netflix. Our
approach supports videos of any resolution, with noadditional latency introduced at the decoder while
supporting random seek.
In addition to the above, we also explored a second
approach which duplicated (rather than moved) the key

3. Proposed Method
Figure 1 illustrates the hierarchical B coding structure
in VVC RA mode where the GOP size is 16 and the Intra
period during encoding is set to 32. Frame-0, which is
encoded as an I-frame, is referenced by multiple frames
within the Intra period. The I-frame is directly referenced
by frames 16, 8, 4, 2, 1 and secondarily referenced byframes 3, 5, 6, 7, 9, 10, 12 and 16 in Picture Order Count(POC). Since these frames are further used in encoding theremaining frames in the same Intra period, the I-frame (i.e.,frame-0) is generally considered to be the most influential
frame and is typically encoded with the aim of minimisingquantisation losses and preserving picture details5).
Our primary method, illustrated in Fig. 2, is based on afull search within the first GOP of the sequence for
optimising RA mode of VVC. The illustration shows anexample of frame-6 as a key frame and is moved to the Iframe position. The new key-frame sequence is thenencoded. Thus, as shown in the illustration, in a GOP sizeof 16, we require 16 key frame sequences to carry out a fullsearch. The aim is to find the best possible key reference
frame that can influence coding efficiency of subsequent
frames, thereby achieving better coding efficiency.
We hypothesise that, by choosing the best key reference
frame from across the entire sequence, further gains can beachieved. However, we limit our search to the first GOP of
the sequence so that there is no increase in the latency at the
decoder. At any time, the decoder has a latency of one GOP
to buffer and reorder the frames according to the POC. The

79

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Table 1 BD-RATE values on 128 sequences from BVI-DVC

reference frame. This increases the frame countcompromising overall BD-RATE gains and hence was not
pursued further. A third method we explored involved
swapping the key reference frame with original frame-0.
Although this approach performs better than simple
duplication, it introduces inter-frame discontinuities anddid not provide any gains due to the additional motion
coding bits required. Due to these frame discontinuities, we
also discounted this method.
In the next section, we discuss the experiments carriedout using the primary method.

Average BD-RATE (%)
Class
Y

U

V

Class A [3940×2160]-3.2522-4.3016-4.832

Class B [1920×1080]-4.2633-5.4201-5.6713Class C [960×540]-4.5982-5.9259-5.7711Class D [480×270]-4.5883-5.7144-5.3904

Overall-4.1755-5.3405-5.4162

4. Experiment Results
BD RATE Y(%)

Our primary approach was validated using the VTM
VVC Software15) version 7.1 in RA mode without any
changes to the codec. We used the default VTM RA
configuration file, again without any changes for either theanchor or our approach. The default RA configuration
disables the rate control function.
We tested our approach on sequences with various
resolutions (the method can be applied to any video
encoded in RA mode). Furthermore, our experiments wereconducted based on the JVET Common Test Conditions7).
The JVET Common Test Conditions lists 26 sequenceswith different classes. We used 13 JVET sequences
corresponding to Class B (HD format), Class C (WVGA
format), and Class D (WQVGA format). We also employed
BVI-DVC sequences to validate results across a wide rangeof video content. The JVET Common Test Conditions
specifies that each sequence must be evaluated against 4
different quantisation values i.e., QP22, QP27, QP32 and
QP37. These values define the quantisation level for the key
frame/I-frame in the sequence. Adaptation and refinementof QP values for the other frames in the GOP are carried out
as specified in VTM RA configuration file.
We evaluated the BD-RATE for each sequence for eachkey frame (16 combinations as described in the previous
section). Additionally, since the luma channel (Y)contributes significantly to the overall picture information,
we calculate BD-RATE on this channel5). The BD-RATE
savings on the chroma channels (U & V) are shown where
possible (cf. Table 1).
We carried out our experiments on a shared cluster witheach node comprising a dual Intel Xeon E5-2670 2.6 GHz,
16 cores in total, 64GB RAM and a 1TB SATA disk
running Linux OS. The facility provides capability torequest nodes on demand specifying number of CPUs andmemory required per node. We conducted our experimentsin parallel processing mode requesting 32GB RAM, 1 CPU
node for each sequence.

2
1.5
1
0.5
0
-0.5

15

14

13

12

11

10

09

08

07

06

05

04

03

02

01

00

Key Frame Selected

BD RATE Y(%)(a) Class B - BasketballDrive

1
0.5
0
-0.5
-1

15

14

13

12

11

10

09

08

07

06

05

04

03

02

01

00

Key Frame Selected(b) Class B - RitualDance

Fig. 3 BD-RATE of Luma channel corresponding to different key
reference frames (a) Class B BasketballDrive (b) Class B
RitualDance

Bjøntegaard Delta (BD-RATE) values are shown in
Table 1. These demonstrate that the technique was able toachieve an average gain of 4.18% (Y channel), which issignificant.
Next, we employed JVET sequences from Class B[1920 ×1080], Class C [832 × 480] and Class D [416 × 240]
sequences for further evaluation of our technique using
encoder intra period size of 64. It should be noted thatsequences in Classes B and C have different resolutions tothose in BVI-DVC. The sequence length used was 81
frames (64 intra period + 1 extra GOP from next intra
period). This enabled us to evaluate the approach acrossIntra periods and the impact of intra period on gains.
An illustration of performance for the JVET Class B
BasketballDrive and RitualDance sequences are shown in
Fig. 3. Different frames were moved and the associated
BD-RATE for the Y-channel is shown for each frame
moved to key frame position. It can be observed that for
BasketballDrive sequence Key Frame 1 gives us the best
spatiotemporal information. Similarly, for RitualDancesequence Key Frame 14 gives us the best spatiotemporal
information with a BD-RATE (Y) saving of 0.75%.
Our experiments indicate that JVET sequences, shownin Table 2, Class C Basketball and PartyScene sequences

4.1 Optimal Key Frame selection on fixed length
sequences
Initially, the frame move was performed on sequencesfrom the BVI-DVC database6) which contains 800 video
sequences with 64 frames each. We randomly selected 128
sequences from this dataset to represent 32 Class A[3940×2160], 32 Class B [1920×1080], 32 Class C[960×540] and 32 Class D sequences [480×270].

80

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Table 2 BD-RATE values on 13 JVET sequences with 81 frames ineach sequence
Class / Sequence[Class B] BasketballDrive

Optimal Key
Frame

BD-RATE (Y)(%)

1-0.206[Class B] BQTerrace

12-0.477[Class B] Cactus

12-0.934[Class B] MarketPlace

15-4.548[Class B] RitualDance14-0.750[Class C] BasketballDrill

0

0.000

[Class C] BQMall

4-0.589[Class C] PartyScene

0

0.000

[Class C] RaceHorses

4-0.921[Class D] BasketballPass

15-2.855[Class D] BlowingBubbles

6-1.956[Class D] BQSquare

15-0.679[Class D] RaceHorses
Average

15

Table 3 BD-RATE values on 13 JVET sequences using frame 16 as
the key frame (81 frames in each sequence)
Class / Sequence-1.395-1.178

BD-RATE (Y) (%)[Class B] BasketballDrive

1.235[Class B] BQTerrace

1.319[Class B] Cactus

0.442[Class B] MarketPlace-5.272[Class B] RitualDance-0.530[Class C] BasketballDrill

4.135[Class C] BQMall

0.908[Class C] PartyScene

3.577[Class C] RaceHorses-0.134[Class D] BasketballPass-2.315[Class D] BlowingBubbles-0.882[Class D] BQSquare-0.722[Class D] RaceHorses-1.086

Average

0.052

as residuals and do not incur large additional bit counts.
Hence, we have been able to save on bits required for
encoding without impacting the residuals of subsequent
frames.
Similarly, Fig. 5 shows the bits used for the Class B
RitualDance sequence with key frames as frame 1, shownin Fig. 5(a), and frame 14, shown in Fig. 5(b). In this
scenario, the bits required for encoding the I-Frame issimilar i.e., approximately 1.2 million bits. However, thebits required for encoding the residuals for POC 1, 2, 3, 13,
16, 32, 48 with key frame 14 are significantly lower whencompared to key frame 1. This confirms our approach, as
described in Section 3, for selecting the best spatiotemporalframe as the I-Frame.

did not benefit from key reference frame selection andhence, in these cases, the anchor is treated as the best-casescenario. Class B sequences showed a BD-RATE (Y)average saving of 1.38%, Class C with 0.38% and Class Dwith 1.72%. The overall saving of 1.2% was achieved. The‘Optimal Key Frame’ for each of JVET sequences is also
shown in Table 2. This demonstrates that the location of thebest possible reference frame is dependent on the sequencecharacteristics i.e., motion involved, motion type (linear or
rotational or combination thereof), amount of backgroundconsistency across frames.
Further, we analysed the bits required for each POC
compared to that for the optimum key frame. With the
optimum key frame, our analysis shows that the bits
required is significantly lower when compared to othercandidate frames. An example is shown in Fig. 4 for theClass B BasketballDrive sequence. Fig. 4(a) shows a casefor the Key Frame 1 while Fig. 4(b) shows a case for the
Key Frame 2. In the charts, the bits required for each POC
are stacked across the QP values i.e., QP 22 is shown inblue, QP 27 is shown in red, QP 32 in green and QP 37 in
purple. It is a known fact QP value specified affects thequantisation process and is a feature that has been aroundfor many years. This means that a lower QP value e.g., QP
22 results in higher picture quality and thus requires
significantly larger number of bits to encode both I-frames
and subsequent B-frames. Similarly, a higher QP value e.g.,
QP 37 produces a lower picture quality and thus requires
lower number of bits. QP 27 and QP 32 are the middle
grounds. This can also be noticed in Fig. 4. We can thus
observe the overall number of bits required to encode each
POC and hence compare how the candidate key frame has
contributed to improving the overall sequence.
From Fig. 4, it can be observed that optimal key frame
1, shown in Fig. 4(a), requires approximately 3 million bits
whereas with key frame 2, Fig. 4(b), approximately 4
million bits are needed. The remaining frames are encoded

4.2 Key Frames at a GOP boundary
We analysed the scenario considering the boundary
condition of a GOP whilst respecting the constraint ofcoding latency. It is well known that VVC (and HEVC)
introduce a latency of 16 frames in RA mode5). The first
frame of the sequence is encoded as an I-frame. The next
frame encoded is frame-16 (GOP size of 16). This
inherently introduces a latency of 17 frames (though frame16 can be treated as belonging to second GOP).
For this scenario, we consider first GOP+1 frame i.e.,
fixing frame 16 as Key Frame. We find that the exceeding
the first GOP of the intra period does not result insignificant savings as shown in Table 3. The only cases
which produce gains for JVET sequences are Class D
BQSquare from –0.679 % to –0.722% in BD-RATE (Y)and Class B MarketPlace from –4.548% to –5.272% in BDRATE (Y).
Comparing results of frame-16 as key frame in Table 3
to the optimal key frame results shown in Table 2, weobserve that all savings are lost. This is due to increase intemporal distance between GOP 0 and GOP 1, seen in Fig.
1. Frame 16 is used as reference in both GOP 0 and GOP 1

81Encoded bits

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Encoded Bits(a) Key Frame Selected = 1(b) Key Frame Selected = 2

Encoded bits

Fig. 4 Comparison of encoded bits for each POC with (a) key frame 1 (optimal key frame) and(b) frame 2 for Class B Basketball sequence

Encoded Bits(a) Key Frame Selected = 1(b) Key Frame Selected = 14

Fig. 5 Comparison of encoded bits for each POC with (a) key frame 1 and(b) frame 14 (optimal key frame) for Class B RitualDance sequenceand thus choosing frame-16 as key frame impacts all Bframes in the intra period.

cSAD(k) is the cumulative SAD for key frame k, cSSD(k)is the cumulative SSD for key frame k and Y[p]x is the pixel
value at position p for frame x.

4.3 Cumulative SAD / SSD values for key frames
We further investigated whether it is possible to reducesearch complexity by determining the optimum key frame
using some metric. Sum of Absolute Difference (SAD) or
Sum of Square Differences (SSD) are typically used as a
metric in image comparison. However, in our scenario, weneed to select the key frame based on spatiotemporalfeatures across the intra period and not just a single frame.
In order achieve spatiotemporal SAD and SSD assessment
across the intra period, we employed a cumulative SAD and
SSD computation. The cumulative SAD and SSD for the Ychannel are described in equations (1) and (2), wherecSAD(k)ൌ  ෍ ෍ ȁܻሾ௣ሿ௙ െ ܻሾ௣ሿ௞ ȁ(1)௙Ǩୀ௞ ௣

cSSD(k)ൌ  ෍ ෍ሺܻሾ௣ሿ௙ െ ܻሾ௣ሿ௞ ሻଶ(2)௙Ǩୀ௞ ௣

The results of the cSAD and cSSD computations areshown in Fig. 6. We mapped the cSAD and cSSD values
against the BD-RATE values observed for the Class B
BasketballDrive (Fig. 6(a)) and RitualDance (Fig. 6(b))
sequences. It can be observed that the optimal key frame for
RitualDance, Fig. 6(b), has similar cSAD and cSSD minima.

82

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Table 4 BD-RATE values on 13 JVET sequences with different sequence lengths
Sequence Length

65 Frames

Class / Sequence

IP Size

Same as IP Size

[Class B] BasketballDrive48

1-0.3481-0.371

1-0.494[Class B] BQTerrace

6412-0.52512-0.82312-0.823[Class B] Cactus

48

9-0.947

9-0.9758-1.509[Class B] MarketPlace

64

15-5.53215-6.67615-6.676[Class B] RitualDance64

14-0.89814-1.07014-1.070[Class C] BasketballDrill

48

0

0.000

0

0.000

0

0.000

[Class C] BQMall

64

4-0.705

5-0.9425-0.942[Class C] PartyScene

64

0

0.000

0

0.000

0

0.000

[Class C] RaceHorses

32

5-1.1895-1.2655-2.555[Class D] BasketballPass

48

15-2.90615-4.02215-4.811[Class D] BlowingBubbles

48

6-2.3366-2.5316-3.233[Class D] BQSquare64

15-0.70715-1.41215-1.412[Class D] RaceHorses

32

15-1.78315-1.7724-3.631

Average BD-RATE (Y) (%)-1.375-1.681-2.089Average BD-RATE (Y) (%)(excluding Class C BasketballDrill and PartyScene)-1.625-1.987-2.469for the RitualDance sequence. Key frame 15 offers thecSAD and cSSD minima, but the BD-RATE shows nobenefit. Hence, we conclude that using cSAD and cSSD to
identify key frames introduces inconsistencies which
deteriorate performance compared to full search. The
investigation of suitable optimisation metrics remains as a
topic for future work.

1

BD RATE Y(%)

2

BD RATE Y(%)

64 Frames

Optimal Key BD-RATE Optimal Key BD-RATE Optimal Key BD-RATE
Frame(Y) (%)Frame(Y) (%)Frame(Y) (%)

1.5
1
0.5
0
-0.5

0.5

0

-0.5

cSAD

4.4 Changing sequence lengths
In the next set of experiments, we set the intra period as
per the JVET test conditions (32 for 30 FPS sequences, 64
for 50 FPS and 60 FPS sequences) with a GOP size of 16.
We selected different sequence lengths as shown in Table
4. For a sequence length of 65 frames, BD-RATE (Y)average savings of 1.4% were obtained, whereas for 64
frames the saving was 1.7%. We can observe that, when thesequence length is the same as the intra period (i.e., 64
frames), additional savings can be achieved. Consequently,
we considered the scenario where the sequence length isidentical to the intra-period length for each sequence (basedon the frame rate). For the case of identical sequence lengthand intra-period, we achieved the best-case BD-RATE (Y)improvement of 2.1%, representing a significant saving
over the conventional VVC encoder.

cSSD

cSSD

cSAD

-1

Key Frame Selected(a) Class B - BasketballDrive

Key Frame Selected(b) Class B - RitualDance

Fig.6 Comparison on BD-RATE seen on Y channel for each key
frame against cSAD and cSSD (Sequence length 81 frames).
(a) Class B BasketballDrive (b) Class B RitualDance

4.5 Workflow, processing times and complexity

However, for the same sequence, we can also see that for
key frame 13, the cSAD and cSSD values though in asimilar region, the BD-RATE is significantly lower.
Figure 6(a) for BasketballDrive, shows an even larger
variation with the cSAD and cSSD values mapped to BDRATE. Only key frame 1 shows a BD-RATE savings while
all others deteriorate for this sequence. Considering thesecSAD and cSSD values, this does not follow the trend seen

We have demonstrated the potential of selecting optimal
key frames using the first GOP of an intra period. The
search evaluates 16 combinations (assuming full search) in
the first GOP. For non-real-time delivery, this search would
be performed in a pre-processing stage (typically in parallelin the cloud) prior to delivery, as shown in Fig. 7(a). Suchparallel pre-processing techniques are widely adopted by83

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）

Create Key
framessequences
Original raw
YUVsequence

Pre-processing: arrange
sequence with X as key frame

sequence with
X = 00 as key frame

Encodersequence with
X = 01 as key frame

Encoder

non-best
bitstreams
bitstreams withbest BD-RATE

sequence with
X = 14 as key frame

Encodersequence with
X = 15 as key frame

Encoder

Compute BD-RATE;Store encoded
bitstreams with best
result; discard rest
selected key frame
metadata

Standard VVC Encoder

Select encoded bitstreams with best BD-Rate(a) Encoder side processing workflow

Fig.7 Workflow and processing times our proposed method(a) Encode side processing workflow (b) Corresponding processing times taken for Class B BasketballDrive

internet streaming providers to prepare a sequence prior to
encoding and delivery.
The 16 key frame sequences are encoded in parallel withno additional processing logic introduced by us at theencoder (or decoder). Once these 16 encodes are complete,
the BD-RATE is calculated for each bitstream. The
bitstream with the best BD-RATE is selected and thenstored together with the key frame position metadata,
shown in Fig. 7(a). The standard compliant bitstream andthe key frame position metadata are used to reconstruct theoriginal picture order at the decoder.
Despite no alterations to the encoder configuration, we
observed a change in the encoder time for different key
frame sequences. As shown in Fig. 7(b), for the Class B
BasketballDrive sequence with frame 00 as key frame, theencoding requires 107,508 seconds, while the sequencewith optimum key frame 01 requires approximately
106,495 seconds. The same sequence with key frame 11
takes 180,272 seconds. This difference is attributed tochange in the spatiotemporal information associated withthe selected key-frame. This spatiotemporal change
influences the motion search across the IntraPeriod (this isa standard VVC algorithm and is not modified by us). It canalso be seen that all other processing time changes are
negligible.
Table 5 shows encoding time required for the optimal
key frame as well as the peak encoding time from the 16
key frame combinations, processed in parallel. It should be
noted that the peak encoding time is equal to the maximum
wait period for the encoder to complete the encoding of 16
key frame sequences.
Our experiment indicates an average overall processing
time overhead at the encoder of approximately 2.2% whenusing the optimal key frame and average of the peakencoder duration of approximately 17. 5%. Further, VVC

Table 5 Increase in processing time at the encoder for optimal key
frame on 13 JVET sequences
Encoder processing
for just the optimal
key frame
Optimal TimeKey
increase
Frame(%)

Class / Sequence

Peak
Encoder
processing
Time
increase(%)[Class B] BasketballDrive1-0.94267.682[Class B] BQTerrace

12

2.757

8.988

[Class B] Cactus

8-0.94320.955[Class B] MarketPlace

15

9.580

10.362[Class B] RitualDance14-1.58224.294[Class C] BasketballDrill

0

0.000

10.883[Class C] BQMall

5

2.673

16.336[Class C] PartyScene

0

0.000

7.965

[Class C] RaceHorses

5-4.045

3.387[Class D] BasketballPass

15

3.698

21.121[Class D] BlowingBubbles

6-1.14914.553[Class D] BQSquare

15

1.190

3.087

[Class D] RaceHorses

4

17.74217.7422.229

17.489Average

RA mode is typically used for delivering pre-encoded
videos and hence, we recommend that our technique isemployed for offline, rather than real time video delivery.

5. Conclusion
In this paper, we have presented a method for improving
the rate-quality performance of VVC based on key frame
selection. We select the best key reference frame from the84

IIEEJ Transactions on Image Electronics and Visual Computing Vol.9 No.2 （2021）(Received Aug. 25, 2021)(Revised Oct. 12, 2021)

first GOP in an intra period, such that it best matches thespatiotemporal information of the sequence, while
maintaining coding latency. Our experiments, using RA
mode, indicate average BD-RATE (Y) savings of 2.1%
over the VTM VVC Software version 7.1 as anchor.
Overall, we consider that the achieved BD-RATE savings,
when compared to the encoder complexity increase, clearlyjustifies the use of our approach.
Our future work will focus on replacing the full key
reference frame search with an efficient algorithm toestimate the best key reference frame.

Sanat Nagaraju
He received his B.E. in Electronics & Communicationfrom Bangalore University in 1999 and has over 20
years of experience with various global telco & media
companies delivering large complex projects in Telco
Networks, IT, and cloud computing. He is currently
pursuing Ph.D. in Visual Communications with theBristol Vision Institute at the University of Bristol.

Fan Zhang
He received the B.Sc. and M.Sc. degrees from
Shanghai Jiao Tong University, Shanghai, China, in2005 and 2008 respectively, and the Ph.D. degree fromthe University of Bristol in 2012. He is a Lecturer in
Visual Communications within the Department ofElectrical and Electronic Engineering, University ofBristol. His research interests include deep videocoding, video quality assessment, perceptual video
compression and creative technology.

Acknowledgement
This research was supported by the NTT Corporation and the
UK Engineering and Physical Sciences Research Council(EPSRC) grant EP/L016656/1 (University of Bristol Centre forDoctoral Training in Communications). This work was carried out
using the computational facilities of the Advanced Computing
Research Centre, University of Bristol –http://www.bristol.ac.uk/acrc/

Seishi Takamura

References
1)

Cisco,
VNI
Complete
Forecast
Highlights,
https://www.cisco.com/c/dam/m/en_us/solutions/serviceprovider/vni-forecasthighlights/pdf/Global_2021_Forecast_Highlights.pdf, p.3, (2021).
2) B. Bross, J. Chen, J. -R. Ohm, G. J. Sullivan, Y. -K. Wang:“Developments in International Video Coding Standardization After
AVC, With an Overview of Versatile Video Coding (VVC)”, Proc. ofIEEE,
Vol.109,
No.9,
pp.1463-1493,
doi:
10.1109/
JPROC.2020.3043399 (2021).
3) G. J. Sullivan, J.-R. Ohm, W. -J. Han, T. Wiegand: “Overview of theHigh Efficiency Video Coding (HEVC) Standard”, IEEE Trans. on
Circuits and Systems for Video Technology, Vol.22, No.12, pp.16491668 (2012).
4) J. Chen, Y. Ye, S. -H. Kim: “Algorithm Description for Versatile
Video Coding and Test Model 7 (VTM 7),” 16th JVET Meeting,
JVET-P2002-v1 (2019).
5) D. R. Bull, F. Zhang, Intelligent Image and Video Compression:Communicating Pictures, 2nd ed., Academic Press Inc. (2021).
6) F. Zhang, D. Ma, D. R. Bull: “BVI-DVC: A Training Database for
Deep
Video
Compression”,
UniversityofBristol,
https://doi.org/10.5523/bris.3hj4t64fkbrgn2ghwp9en4vhtn (2020).
7) J. Boyce, K. Suehring, X. Li, V. Seregin: “JVET Common Test
Conditions and Software Reference Configurations”, 10th JVET
Meeting, JVET-J1010-v1 (2018).
8) S. De-Luxán-Hernández, G. Valeri, M. Jackie, N. Tung, S. Heiko, M.
Detlev, T. Wiegand: “An Intra Subpartition Coding Mode for VVC”,
Proc. of IEEE International Conference on Image Processing(ICIP2019), pp.1203-1207 (2019).
9) G. Bjøntegaard: “Calculation of Average PSNR Differences Between
RD-Curves”, 13th VCEG Meeting, VCEG-M33 (2001).
10) G. Bjøntegaard: “Improvements of the BD-PSNR Model”, ITU-T
SG16/Q6, 35th VCEG Meeting, VCEG-AI11 (2008).
11) H. Schwarz, D. Marpe, T. Wiegand: “Analysis of Hierarchical B
Pictures and MCTF”, Proc. of IEEE International Conference on
Multimedia
andExpo(ICME2006),
pp.1929-1932,
doi:10.1109/ICME.2006.262934 (2006).
12) B. Erabadda, T. Mallikarachchi, G. Kulupana, A. Fernando: “Virtual
Frames as Long-Term Reference Frames for HEVC Inter-Prediction”,
Proc. of IEEE International Conference on Consumer Electronics(ICCE2020), pp.1-2 (2020).
13) S. Takamura, A. Shimizu: “Simple and Efficient H.265/HEVC
Coding of Fixed Camera Videos”, Proc. of IEEE International
Conference on Image Processing (ICIP2016), pp.804-808 (2016).
14) HEVC Reference Software (HEVC Test Model HM16.4)
https://vcgit.hhi.fraunhofer.de/jvet/HM (2021).
15) VVC Reference Software (VVC Test Model VTM7.1),
https://vcgit.hhi.fraunhofer.de/jvet/VVCSoftware_VTM//tree/VTM-7.1 (2020).

85(Member)
He received the B.E., M.E. and Ph.D. degrees from the
University of Tokyo in 1991, 1993 and 1996,
respectively. In 1996 he joined Nippon Telegraph andTelephone (NTT) Corporation, where he has been
engaged in research on efficient video coding andmulti-modal signal modelling. He has fulfilled various
duties in the research and academic community incurrent and prior roles including Associate Editor ofIEEE Transactions on Circuits and Systems for Video
Technology, Executive Committee Member of IEEE
Tokyo Section, Japan Council and Region 10 (AsiaPacific). He has also served as Chair of ISO/IEC JTC
1/SC 29 Japan National Body, Japan Head ofDelegation of ISO/IEC JTC 1/SC 29, and as an
International Steering Committee Member of the
Picture Coding Symposium. From 2005 to 2006, he
was a Visiting Scientist at Stanford University,
California, USA. He is currently a SeniorDistinguished Engineer at NTT Computer and Data
Science Laboratories. He is a Fellow of IEEE, IEICE
and ITE, a senior member of IPSJ, and a member of
MENSA, APSIPA (VP-IRD), and SID.

David R. Bull
He received the B.Sc. degree from the University of
Exeter, Exeter, U.K., in 1980, the M.Sc. degree fromthe University of Manchester, Manchester, U.K., in
1983, and the Ph.D. degree from the University ofCardiff, Cardiff, U.K., in 1988. He has previously beena Systems Engineer with Rolls Royce, Bristol, U.K.
and a Lecturer with the University of Wales, Cardiff,
U.K. He was with the University of Bristol in 1993 andis currently its Chair of Signal Processing and Director
with Bristol Vision Institute. He is also Director withthe recently announced £46m UKRI ‘MyWorld’Strength in Places Programme. In 2001, he co-foundedauniversity spin-off
company,
ProVisionCommunication Technologies Ltd., specializing inwireless video technology. He has authored or coauthored more than 500 papers on the topics of image
and video communications and analysis for wireless,Internet and broadcast applications, together withnumerous patents, several of which have been
exploited commercially. He is the author of threebooks and has delivered numerous invited/keynote
lectures and tutorials. He was the recipient of the two
IET Premium awards for his work. He is a Fellow ofIEEE and the Institution of Engineering andTechnology.

Call for Papers
6SHFLDO,VVXHRQWKH([WHQGHG3DSHUV3UHVHQWHGLQWKHWK,,((-,QWHUQDWLRQDO
&RQIHUHQFHRQ,PDJH(OHFWURQLFVDQG9LVXDO&RPSXWLQJ ,(9&                                  ,,((-(GLWRULDO&RPPLWWHH
7KHWK,,((-,QWHUQDWLRQDO&RQIHUHQFHRQ,PDJH(OHFWURQLFVDQG9LVXDO&RPSXWLQJ ,(9& ZDVKHOGRQOLQHRQ6HSWHPEHU7KHDLPRIWKHFRQIHUHQFHLVWREULQJWRJHWKHUUHVHDUFKHUVHQJLQHHUVGHYHORSHUV DQG VWXGHQWV IURP YDULRXV ILHOGV LQ ERWK DFDGHPLD DQG LQGXVWU\ IRU GLVFXVVLQJ WKH ODWHVWUHVHDUFKHV VWDQGDUGV GHYHORSPHQWV LPSOHPHQWDWLRQV DQG DSSOLFDWLRQ V\VWHPV LQ DOO DUHDV RI LPDJHHOHFWURQLFVDQGYLVXDOFRPSXWLQJ7KH&RQIHUHQFHDOUHDG\VROLFLWHG-RXUQDO7UDFN -7 3DSHUVWREHLQFOXGHG
LQWKH6SHFLDO,VVXHRQ³-RXUQDO7UDFN3DSHUVLQ,(9&´WREHSXEOLVKHGRQ'HFHPEHULVVXHRIWKH
7UDQVDFWLRQVRQ,PDJH(OHFWURQLFVDQG9LVXDO&RPSXWLQJ 9RO1R  
7KH HGLWRULDO FRPPLWWHH DOVR SODQV WR SXEOLVK WKH VSHFLDO LVVXH RQ ³([WHQGHG 3DSHUV 3UHVHQWHG LQ
,(9&´ RQ -XQH  LVVXH RI WKH ,,((- 7UDQVDFWLRQV RQ ,PDJH (OHFWURQLFV DQG 9LVXDO &RPSXWLQJ
9RO1R 7KHHGLWRULDOFRPPLWWHHZLOOZLGHO\DVNIRUVXEPLVVLRQVIRUWKHSDSHUVRIWKHIROORZLQJDUHD
E\ H[WHQGLQJ WKH PDWHULDO SUHVHQWHG LQ DOO VHVVLRQV RI ,(9& :H KRSH \RX ZRXOG VXEPLW \RXUKLJKTXDOLW\RULJLQDOSDSHUVDIWHUFKHFNLQJWKHHOHFWURQLFVXEPLVVLRQJXLGHOLQHVRQRXUVLWH 7RSLFVFRYHUHGLQFOXGHEXWDUHQRWOLPLWHGWR,PDJH DQG 9LGHR &RGLQJ 7UDQVFRGLQJ &RGLQJ 6WDQGDUGV DQG 5HODWHG ,VVXHV ,PDJH DQG 9LGHR
3URFHVVLQJ ,PDJH $QDO\VLV 6HJPHQWDWLRQ DQG &ODVVLILFDWLRQ ,PDJH 5HFRJQLWLRQ ,PDJH 5HVWRUDWLRQ
6XSHU5HVROXWLRQ &RORU 5HVWRUDWLRQ DQG 0DQDJHPHQW &RPSXWHU 9LVLRQ 0RWLRQ $QDO\VLV &RPSXWHU*UDSKLFV 0RGHOLQJ 5HQGHULQJ 9LVXDOL]DWLRQ $QLPDWLRQV ,QWHUDFWLRQ 135 9LUWXDO 5HDOLW\ DQG ',PDJLQJ 'DWD +LGLQJ :DWHUPDUNLQJ DQG 6WHJDQRJUDSK\ &RQWHQW 3URWHFWLRQ %LRLQIRUPDWLFV DQG$XWKHQWLFDWLRQ&RPSXWHU)RUHQVLFV,PDJHGDWDEDVH,PDJHDQG9LGHR5HWULHYDO'LJLWDOPXVHXP'LJLWDO$UFKLYLQJ &RQWHQW 'HOLYHU\ ,PDJH $VVHVVPHQW ,PDJH 4XDOLW\ 3ULQWLQJ DQG 'LVSOD\ 7HFKQRORJLHV,PDJLQJ 'HYLFHV 'LJLWDO 6LJQDJH (OHFWURQLF 3DSHU 9LVXDO &RPPXQLFDWLRQ +XPDQ ,QWHUIDFHV DQG,QWHUDFWLRQV 0RELOH ,PDJH &RPPXQLFDWLRQ 1HWZRUNLQJ DQG 3URWRFROV 2SWLFDO &RPPXQLFDWLRQ+DUGZDUH DQG 6RIWZDUH ,PSOHPHQWDWLRQ ,PDJH 5HODWHG $SSOLFDWLRQV /6, 8QGHUVWDQGLQJ RI +XPDQ
9LVLRQ DQGRU +XPDQ 7DFWLOH 6HQVH :HE5HODWHG 7HFKQLTXHV 3HUVRQDOL]DWLRQ 7HFKQLTXH ,QWHUDFWLRQ
EHWZHHQ+XPDQDQG&RPSXWHU8VDELOLW\$FFHVVLELOLW\,PDJH3URFHVVLQJ7HFKQLTXH&RQVLGHUHG(PRWLRQ2WKHU)XQGDPHQWDODQG$SSOLFDWLRQ7HFKQLTXH,QWHUQDWLRQDO6WDQGDUGL]DWLRQ 7UHDWPHQWRISDSHUV6XEPLVVLRQSDSHUVW\OHIRUPDWDQGGRXEOHEOLQGSHHUUHYLHZSURFHVVDUHWKHVDPHDVWKHUHJXODUSDSHU 
 3XEOLFDWLRQRI6SHFLDO,VVXH,((-7UDQVDFWLRQVRQ,PDJH(OHFWURQLFVDQG9LVXDO&RPSXWLQJ9RO1R -XQH  
 6XEPLVVLRQ'HDGOLQH
)ULGD\)HEUXDU\
 &RQWDFWGHWDLOVIRU,QTXLUHV ,,((-2IILFH(PDLOKHQV\X#LLHHMRUJ
 2QOLQH6XEPLVVLRQ85/KWWSZZZHGLWRULDOPDQDJHUFRPLLHHM 

86

Call for Papers
6SHFLDO,VVXHRQ'DWD,QWHUZRUN6KDULQJDQG5HXVLQJ7HFKQRORJLHV6XSSRUWLQJ'LJLWDO7UDQVIRUPDWLRQ(UD,,((-(GLWRULDO&RPPLWWHH
,QRUGHUWRVROYHYDULRXVVRFLDOLVVXHVVSHHG\DQGHFRQRPLFDOO\'LJLWDO7UDQVIRUPDWLRQ '; LVH[SHFWHG
WR UHDOL]H )XOO\ 'LJLWDOL]HG 6RFLHW\ 7KHUHIRUH QRW RQO\ WKH QDWXUDO ODQJXDJH DQG LPDJHYLGHR RI WKH RXWHUZRUOGLQGLYLGXDOEHKDYLRUDQGHYHQWKHPRWLYDWLRQDQG.DQVHLRIWKHVROHKXPDQDUHJRLQJWREHSURFHVVHGGLJLWDOO\7KHFRQFHSWRIVPDUWKRPHDQGVPDUWFLW\ZLOOEHWKHH[DPSOHRIVXFKWUHQG7KHXWLOL]DWLRQRIWKH
GDWDDXWRPDWLFDOO\DFTXLUHGIURPWKHXVHUZLOOLPSURYHWKHVSHFLILFDWLRQRIWKHSURGXFWVDQGVKDULQJVHUYLFHZLOO GHJUDGH WKH KXUGOH RI SRVVHVVLRQ RI WKH H[SHQVLYH SURGXFWV IRU WKH FXVWRPHUV LQ VXFK SHUVSHFWLYH 7R
SURPRWHVXFKWLGHWKHLQWHUZRUNVKDULQJDQGUHXVLQJWKHGLJLWDOGDWDZLWKRXWFRQYHUWLQJLQWRDQDORJXHGDWDRUZLWKRXWWUXQFDWLQJLWLVYHU\LPSRUWDQW
(VSHFLDOO\LQLPDJHUHODWHGDUHDWKHWHFKQRORJLHVWRFRQYHUWDQGSURFHVVWKHLPDJHGDWDDQGLWVDWWULEXWHVZLOOFRQWULEXWHWRSURGXFHWKHYLUWXDORXWORRNRIWKHSURGXFWVZLWKRXWPDQXIDFWXULQJDQGWKHLPDJHDQDO\VLV
E\$,ZLOOSURYLGHIXUWKHUEHQHILWLQEXVLQHVV
7KLV VSHFLDO LVVXH RI WKH SDSHU WDUJHWV YDULRXV LPDJHUHODWHG WHFKQRORJLHV VXFK DV LPDJH SURFHVVLQJ
JHQHUDWLRQ H[SUHVVLRQ WUDQVPLVVLRQ DQG FRPELQHG DSSOLFDWLRQV RI WKHVH WKDW VXSSRUW DQG SURPRWH GDWD
LQWHUZRUNVKDULQJDQGUHXVLQJIRUWKHGLJLWDOWUDQVIRUPDWLRQHUDDQGFDOOVIRUDQ\FDWHJRU\ 2UGLQDU\SDSHU6KRUWSDSHU6\VWHPGHYHORSPHQWSDSHU'DWDSDSHU3UDFWLFH2ULHQWHG3DSHU RISDSHUV  7RSLFVFRYHUHGLQFOXGHEXWQRWOLPLWHGWR ,PDJH 3URFHVVLQJ ,PDJH 5HFRJQLWLRQ ,PDJH 'HWHFWLRQ ,PDJH9LGHR &RPPXQLFDWLRQ ,PDJH,QSXW2XWSXW &RPSXWHU *UDSKLFV 9LVXDOL]DWLRQ %LQRFXODU 9LVLRQ ' LPDJH SURFHVVLQJ &RPSXWHU9LVLRQ 3DWWHUQ 5HFRJQLWLRQ %LJ 'DWD ,PDJH 'DWDEDVHV 0DFKLQH /HDUQLQJ 'HHS /HDUQLQJ
8QGHUVWDQGDELOLW\  $QQRWDWLRQ $WWULEXWH 3UHVXPSWLRQ :HEEDVHG 7HFKQRORJ\ 6HFXULW\ &UHDWLYLW\8VDELOLW\ ,QWHUSUHWDELOLW\ +XPDQ ,QWHUIDFH DQG ,QWHUDFWLRQ 8VHU ([SHULHQFH 8ELTXLWRXV $59505
2WKHUUHODWHGIXQGDPHQWDODSSOLFDWLRQV\VWHPL]HGWHFKQRORJLHV 7UHDWPHQWRISDSHUV6XEPLVVLRQSDSHUVW\OHIRUPDWDQGGRXEOHEOLQGSHHUUHYLHZSURFHVVDUHWKHVDPHDVWKHUHJXODUSDSHU,I
WKH QXPEHU RI DFFHSWHG SDSHUV LV OHVV WKDQ WKH PLQLPXP QXPEHU IRU WKH VSHFLDO LVVXH WKH DFFHSWDQFH
SDSHUZLOOEHSXEOLVKHGDVWKHUHJXODUFRQWULEXWHGSDSHU:HDVNIRU\RXUXQGHUVWDQGLQJDQGFRRSHUDWLRQ 
 3XEOLFDWLRQRI6SHFLDO,VVXH,((-7UDQVDFWLRQVRQ,PDJH(OHFWURQLFVDQG9LVXDO&RPSXWLQJ9RO1R 'HFHPEHU  
 6XEPLVVLRQ'HDGOLQH
)ULGD\-XO\
 &RQWDFWGHWDLOVIRU,QTXLUHV ,,((-2IILFH(PDLOKHQV\X#LLHHMRUJ
 2QOLQH6XEPLVVLRQ85/KWWSZZZHGLWRULDOPDQDJHUFRPLLHHM 

87

Call for Papers
6SHFLDO,VVXHRQ,PDJH5HODWHG7HFKQRORJ\IRU'LJLWDO)DEULFDWLRQ,,((- (GLWRULDO &RPPLWWHH,QUHFHQW\HDUVWKHUHKDVEHHQDORWRIUHVHDUFKRQGLJLWDOIDEULFDWLRQZKLFKLVWKHSURFHVVRI
FUHDWLQJVRPHWKLQJEDVHGRQGLJLWDOGDWDDQGYDULRXVIDEULFDWLRQWRROVVXFKDV'SULQWHUVODVHUFXWWHUV'VFDQQHUV&1&PLOOLQJPDFKLQHVHWFKDYHEHFRPHDFFHVVLEOHWRWKHJHQHUDOSXEOLF
HQDEOLQJ LQGLYLGXDOV WR FUHDWH QHZ WKLQJV $V D UHVXOW LW KDV EHFRPH SRVVLEOH WR FUHDWH QHZSURGXFWVRQDSHUVRQDOOHYHO,QDGGLWLRQWRSODVWLFPRGHOLQJPHWDODQGZRRGDUHDOVRDYDLODEOH
DVPDWHULDOVDQGWKHXVHRIYDULRXVPDWHULDOVLQFRPELQDWLRQVXFKDVSULQWLQJHOHFWURQLFFLUFXLWV
KDV EHHQ SURSRVHG DQG WKH GHYHORSPHQW DQG DSSOLFDWLRQ RI 0(06 LV DOVR H[SHFWHG 7KXVGLJLWDOIDEULFDWLRQLVQRWRQO\DUDSLGSURWRW\SLQJRISURGXFWVEXWDOVRDWHFKQRORJ\WKDWFDQEHDGRSWHGLQDZLGHUDQJHRIILHOGVVXFKDVIRRGFORWKLQJDUFKLWHFWXUHPHGLFLQHHGXFDWLRQDQG
HQWHUWDLQPHQW
,Q WKLV VSHFLDO LVVXH ZH LQYLWH FDWHJRU\ 2UGLQDU\ SDSHU 6KRUW SDSHU 6\VWHP GHYHORSPHQWSDSHU'DWDSDSHU3UDFWLFH2ULHQWHG3DSHU RISDSHUVWRUHDOL]HGLJLWDOIDEULFDWLRQ:HZRXOGEH
JUDWHIXOLI\RXFRXOGVXEPLW\RXUSDSHUV7RSLFVFRYHUHGLQFOXGHEXWQRWOLPLWHGWR&RPSXWHU *UDSKLFV ,PDJH 3URFHVVLQJ ,PDJH $QDO\VLV ,PDJH 5HFRJQLWLRQ ,QWHUDFWLRQ&RPSXWHU 9LVLRQ 0DFKLQH /HDUQLQJ'HHS /HDUQLQJ 8VHU ,QWHUIDFH 95 $5 05 '
5HODWHG 7HFKQRORJLHV ,PDJH 6\VWHPV ,PDJH $SSOLFDWLRQV *HQHUDO ,PDJH %DVLF
7HFKQRORJLHV7UHDWPHQWRISDSHUV6XEPLVVLRQ SDSHU VW\OH IRUPDW DQG GRXEOHEOLQG SHHU UHYLHZ SURFHVV DUH WKH VDPH DV WKHUHJXODU SDSHU ,I WKH QXPEHU RI DFFHSWHG SDSHUV LV OHVV WKDQ WKH PLQLPXP QXPEHU IRU WKH
VSHFLDO LVVXH WKH DFFHSWDQFH SDSHU ZLOO EH SXEOLVKHG DV WKH UHJXODU FRQWULEXWHG SDSHU :H
DVNIRU\RXUXQGHUVWDQGLQJDQGFRRSHUDWLRQ 
3XEOLFDWLRQRI6SHFLDO,VVXH,((-7UDQVDFWLRQVRQ,PDJH(OHFWURQLFVDQG9LVXDO&RPSXWLQJ9RO1R -XQH 
6XEPLVVLRQ'HDGOLQH
0RQGD\2FWREHU
&RQWDFWGHWDLOVIRU,QTXLUHV ,,((-2IILFH(PDLOKHQV\X#LLHHMRUJ
2QOLQH6XEPLVVLRQ85/KWWSZZZHGLWRULDOPDQDJHUFRPLLHHM 

88

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:z 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.
z 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.
z 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.
z 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 Word89

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:
9 Novelty: The contents of the paper are novel.
9 Utility: The contents are useful for academic and industrial development.
9 Reliability: The contents are considered trustworthy by the reviewer.
9 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.
9 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.
9 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 (as90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.

91

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

Vol. 9, No. 2 2021

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