Takahiko Furuya and Ryutarou Ohbuchi, Diffusion-on-Manifold Aggregation of Local Features for Shape-based 3D Model Retrieval,
accepted as oral paper, Proc. ACM International Conference on Multimedia Information Retrieval (ICMR) 2015, Shanghai, China.

Diffusion-on-Manifold Aggregation of Local Features
for Shape-based 3D Model Retrieval
Takahiko Furuya

Ryutarou Ohbuchi

University of Yamanashi
4-3-11 Takeda, Kofu-shi
Yamanashi-ken, 400-8511, Japan+81-55-220-8570

University of Yamanashi
4-3-11 Takeda, Kofu-shi
Yamanashi-ken, 400-8511, Japan+81-55-220-8570

g13dm003AT yamanashi.ac.jp

ohbuchiAT yamanashi.ac.jp

ABSTRACT

1. INTRODUCTION

Aggregating a set of local features has become one of the mostcommon approaches for representing a multi-media data such as
2D image and 3D model. The success of Bag-of-Features (BF)
aggregation [2] prompted several extensions to BF, that are, VLAD[12], Fisher Vector (FV) coding [22] and Super Vector (SV) coding[34]. They all learn small number of codewords, or representative
local features, by clustering a set of large number of local features.
The set of local features extracted from a media data (e.g., animage) is encoded by considering distribution of features aroundthe codewords; BF uses frequency, VLAD and FV uses
displacement vector, and SV uses a combination of both. In doing
so, these encoding algorithms assume linearity of feature spaceabout a codeword. Consequently, even if the set of features form a
non-linear manifold, its non-linearity would be ignored, potentiallydegrading quality of aggregated features. In this paper, we propose
a novel feature aggregation algorithm called Diffusion-on-Manifold(DM) that tries to take into account, via diffusion distance, structureof non-linear manifold formed by the set of local features. In viewof 3D shape retrieval, we also propose a local 3D shape featuredefined for oriented point set. Experiments using shape-based 3Dmodel retrieval scenario show that the DM aggregation results inbetter retrieval accuracy than the existing aggregation algorithms
we’ve compared against, that are, VLAD, FV, and SV, etc.

Majority of state-of-the-art shape-based 3D model retrieval(3DMR) algorithms represent a 3D model as a set of local features
extracted from its 3D geometry or 2D appearances [4, 5, 6, 14, 15,
17, 19, 20]. Using a set of local features with an appropriate
aggregation algorithm yields a 3D shape descriptor robust against
articulation or global deformation of 3D models. Aggregation of
local feature also make comparison among a pair of 3D models
much more efficient than comparing local features individually.
Bag-of-Features (BF) [2] (Figure 1a) is first-of-a-kind local feature
aggregation algorithm and is widely used to realize classification,
segmentation, annotation, or retrieval of 2D images or 3D models.
BF first learns a set of codewords by clustering the set of local
features extracted from all the 3D models in a training set database.
Each of a set of local features extracted from a 3D model is vectorquantized into its nearest codeword, and the codewords areaccumulated into a histogram to become an aggregated feature
vector for the 3D model. The vector quantization (VQ) process of
the BF throws away a lot of information, e.g., relative positions or
global arrangement of the set of local features. Consequently, thequality of aggregated feature suffers.
Several refinements to BF aggregation have been proposed [12, 16,
22, 23, 29, 30, 31, 34] to remedy this issue. They use either sparse
coding (SC) or higher-order statistics (HS) to reduce information
loss due to VQ. Both SC-based and HS-based algorithms use anumber of codewords much smaller than the one for BF, andcomplement VQ with additional information. The Local
Coordinate (LC) coding [31], Locality-constrained Linear (LL)
coding [29], and localized soft-assigment coding [16] sparsely
encode each local feature as weighted linear sum of neighboring
codewords (Figure 1b). The Vector of Locally Aggregated
Descriptor (VLAD) [12], Fisher Vector (FV) coding [22], andSuper Vector (SV) coding [34] all try to encode local features byusing higher-order statistics around the codewords (Figure 1c). The
FV aggregates the set of local features by using displacement ofmean and variance between the local features and the codewords,
each of which is defined as a multivariate normal distribution
derived by GMM clustering. The VLAD, which can be considered
as an approximation of the FV, uses displacement of mean only.

Categories and Subject Descriptors
H.3.3 [Information Search and Retrieval]: Information filtering.
I.2.10 [Vision and Scene Understanding]: 3D/stereo scene
analysis.

General Terms
Algorithms, Experimentation.

Keywords
Feature aggregation, content-based retrieval, manifold learning,
manifold ranking, diffusion distance, local feature, bag-of-features,sparse coding, fisher vector, super vector, 3D shape, 3D oriented
point set.
Permission to make digital or hard copies of all or part of this work forpersonal or classroom use is granted without fee provided that copies are
not made or distributed for profit or commercial advantage and that copies
bear this notice and the full citation on the first page. Copyrights forcomponents of this work owned by others than ACM must be honored.
Abstracting with credit is permitted. To copy otherwise, or republish, topost on servers or to redistribute to lists, requires prior specific permission
and/or a fee. Request permissions from Permissions@acm.org.

These SC-based and HS-based algorithms aggregate a set of local
features more accurately than a VQ-based approach, i.e., BF.
However, most of them still can’t handle non-linear distribution of
local features in feature space. Small number of codewords
essentially ignores local geometry of feature distribution.
Extensions employed by the SC-based and HS-based algorithms
listed above won’t be able to completely capture non-linear natureof distribution of a set of local features in its feature space.

ICMR '15, June 23 - 26, 2015, Shanghai, China
Copyright 2015 ACM 978-1-4503-3274-3/15/06…$15.00
http://dx.doi.org/10.1145/2671188.2749380

1Takahiko Furuya and Ryutarou Ohbuchi, Diffusion-on-Manifold Aggregation of Local Features for Shape-based 3D Model Retrieval,
accepted as oral paper, Proc. ACM International Conference on Multimedia Information Retrieval (ICMR) 2015, Shanghai, China.
An experimental evaluation of proposed DM aggregation using
multiple local features and multiple benchmark databases showsthat it significantly outperforms existing SC-based and HS-based
aggregation methods. Furthermore, POD features aggregated byusing DM yields retrieval accuracy comparable to state-of-the-art
3DMR algorithms using local 2D visual features.

local feature

codeword(a) VQ-based aggregation.

Contributions of this paper can be summarized as follows.
(b) SC-based aggregation.

 Proposition of Diffusion-on-Manifold (DM) aggregation of
local features. DM exploits non-linear manifold structure of
local features for more accurate aggregation.

manifold graphof local features

 Proposition of a local 3D geometric feature for oriented point
set. POD feature performs comparably in retrieval accuracy
to state-of-the-art local 2D visual features for 3DMR.

(c) HS-based aggregation.

The rest of this paper is structured as follows. We will describe
related work in the next section. The proposed algorithms aredescribed in Section 3. Empirical evaluation of the proposed
algorithms will be presented in Section 4, followed by conclusionand future work in Section 5.

(d) Proposed DM aggregation.

Figure 1. The proposed Diffusion-on-Manifold (DM)algorithm (d) aggregates a set of local features by relevancediffusion on the manifold graph of local features to generate
more accurate aggregated features than VQ-based (a), SCbased (b), and HS-based (c) approaches.

2. RELATED WORK
2.1 Aggregation of Local Features
The BF [2] is first-of-a-kind feature aggregation algorithm and was
first introduced to 3DMR by Liu et al. [17]. The BF counts thenumber of vector-quantized local features to generate frequency
histogram of codewords. Retrieval accuracy of BF-aggregatedfeatures are insufficient since positional information of local
features in their feature space is lost due to VQ.

In this paper, we propose a novel local feature aggregation
algorithm called Diffusion-on-Manifold (DM) that exploits
structure of potentially non-linear manifold of local features. Figure
1 illustrates the proposed algorithm along with the existingalgorithms for feature aggregation. Previous aggregation
algorithms learn a set of cluster centers, that are, codewords, froma large number (e.g., 250k) of local features in a training set. During
training stage, our algorithm generates a graph representing
potentially non-linear manifold of all the training set features.
Given a query, e.g., a 3D model, unseen local features extracted
from it are vector quantized to their respective nearest neighbors in
the manifold graph. Then, using these local features of the query assources, relevance values are diffused over the manifold graph.
After diffusion, relevance values from all the features on the
manifold graph are aggregated into a feature vector for the query.
The proposed DM aggregation algorithm may be regarded as anextension of LL coding [29] to non-linear manifold.

To more accurately aggregate the set of local features, several SCbased aggregation algorithms and HS-based aggregation
algorithms have been proposed. The ScSPM [30], LC coding [31],
LL coding [29], and localized soft-assignment coding [16] aregrouped in SC-based approach. They sparsely encode each local
feature as weighted linear sum of neighboring codewords. The
VLAD [12], VLAT [23], FV [22], and SV [34] are grouped in HSbased approach. They all try to encode local features by using
higher-order statistics around the codewords. They computedisplacement vector between the mean of local features and thecodeword, possibly with some additional information. In addition
to displacement vector of mean, FV accumulates displacement
vector with respect to variance of local features and SV
accumulates frequency of codewords as with the BF. VLAT
encodes local features around the codeword by using theircovariance in addition to displacement vector of mean.

In view of 3DMR, in addition to the feature aggregation algorithm,
the choice of low-level local shape feature to be aggregated isequally important. Local features for 3D models can be classified
into two groups; local 3D geometric feature [8, 13, 17, 20, 24, 27]
and local 2D visual feature [4, 5, 6, 19]. At this point, the latter hasan advantage in retrieval accuracy. Many state-of-the-art 3DMR
algorithms employ local 2D visual features extracted densely fromimages rendered from multiple viewpoints of 3D model. For
example, top two finishers in the SHape REtrieval Contest(SHREC) 2014 Large-scale Comprehensive 3D shape retrieval(SH14LC) track used aggregation of local 2D visual features [15].
Despite their success, a 2D visual feature has its drawback; it isincapable of describing internal structures of 3D models. To
capture internal structure of 3D models represented as surfaces,
voxels, or point set, an accurate and efficient local 3D geometricalfeature is desired.

These SC-based approaches and HS-based approaches produce
higher accuracy than BF in 2D image classification [1] or 3DMR
[5]. However, accuracy of these methods may be still insufficientsince they ignore structure of non-linear manifold of local features.
The FV can essentially capture non-linearity of features due to its
use of multivariate Gaussian distributions as codewords. But, inpractice, it is difficult to approximate the structure of feature
manifold by using small number (e.g., tens to hundreds) ofmultivariate Gaussians.

2.2 Manifold Learning
Manifold learning is used in many applications such as information
retrieval, classification, clustering, or segmentation, to improvetheir accuracies. It learns better distance metric in the feature spaceby analyzing the structure of non-linear distribution of features, i.e.,
manifold. Typically, the manifold is represented as a graph where
each node corresponds to a feature of a multi-media data (e.g., 3D
model), and each edge indicates similarity between two features.

We thus propose a local 3D geometric feature called Position andOrientation Distribution (POD) for 3D oriented point set. To
accurately describe a geometry represented by distribution of a setof oriented points enclosed in a local region of a 3D model, POD
encodes the oriented point sets by using SV coding.
2

Takahiko Furuya and Ryutarou Ohbuchi, Diffusion-on-Manifold Aggregation of Local Features for Shape-based 3D Model Retrieval,
accepted as oral paper, Proc. ACM International Conference on Multimedia Information Retrieval (ICMR) 2015, Shanghai, China.
weights are used in relevance diffusion on the manifold graph formore accurate aggregation.

For information retrieval, several diffusion distance-based ranking
algorithms have been proposed [3] to improve retrieval accuracy.
One of the representative algorithms is the Manifold Ranking (MR)by Zhou et al. [33]. The MR diffuses relevance from source node(s),which correspond to given query(s), over the manifold graph.
Retrieval ranking is generated base on relevance values diffused on
the manifold.

The POD feature describes distribution of oriented points within a
SOI by using SV coding. The SOI is divided into multiple cells byusing a regular grid. For each cell, positions of the oriented points
within the cell are encoded by using SV, and normals of theoriented points are described by their covariance. A set of POD
features is densely extracted from a 3D model, and is aggregatedby the DM algorithm for accurate 3D model comparison.

While most of studies using manifold learning apply it on thefeature space of multi-media data such as 3D model, only a few
studies apply it on the local feature space. Tao et al. [26] computed
the MR on the manifold of local 3D geometric features extracted
from a 3D model for mesh saliency detection. Torki et al. [28] andZhu et al. [35] reduced dimensionality of local features by eigenanalyzing the manifold graph of local 2D visual features for
accurate feature matching. In this paper, we utilize the manifold of
local features for feature aggregation.

3.2 Diffusion-on-Manifold Aggregation
3.2.1 Aggregating Local Features
Given a set of local features extracted from a 3D model, DM
algorithm aggregates the set of local features into a feature vector
per 3D model by relevance diffusion on the manifold graph of local
features. The manifold graph is represented as a sparse matrix P,
whose construction process will be described in Section 3.2.2.1.
The matrix P has a size Nt×Nt where Nt is the number of nodes, i.e.,
training local features, on the manifold graph (e.g., Nt=250K).

2.3 3D Local Feature for Oriented Point Set
Recently, 3D oriented point set has been becoming one of the mostused 3D shape representation due to proliferation of inexpensive
3D range scanners. To compare or recognize 3D models
represented as oriented point set, many local 3D geometric features
for oriented point set have been proposed [8, 13, 20, 24, 27].

The DM algorithm first generates an Nt-dim. source vector y forrelevance diffusion. For each local feature x of a 3D model, x is
vector-quantized into its nearest node (i.e., nearest training local
feature) and the nearest node is set as a source for relevancediffusion. We use a kd-forest having 20 kd-trees to efficiently andaccurately search nearest nodes. The source value yn for the node nin the source vector y is computed by using the following equation;

Given an oriented point set of a 3D model, a set of Sphere-OfInterest (SOI), in which a local feature is computed, is densely
sampled. SPFH [24] describes the SOI by angular statistics of pairsof oriented points to form a 125-dim. feature vector. LSF [20] uses
distances among the oriented points within the SOI in addition toangular statistics to generate 625-dim. feature vector. Spin Image(SI) [13] projects a set of oriented points in the SOI into acylindrical coordinate and counts frequency of points for eachregion of the coordinate. SI feature is a 153-dim. vector. SPFH, LSF,
and SI are invariant against 3D rotation of SOIs. RoPS [8] describes
density of the point set within the SOI. It first normalizes 3D
rotation of the SOI by using PCA. The rotation-normalized SOI isthen rendered from multiple viewpoints to generate a set of 2D
images. For each rendered image, RoPS computes first- andsecond-order moments to generate 135-dim. feature vector.

yn  wn N n(1)

where Nn is the number of local features vector-quantized into thenode n, and wn is the weight for the node n. An algorithm forcomputing wn will be described in Section 3.2.2.2.
We then compute relevance diffusion on the manifold graph togenerate the aggregated feature vector for the 3D model. Relevanceis diffused from the multiple sources defined by the source vector
y over the manifold graph P. We use the following iterative formof relevance diffusion used in the Manifold Ranking algorithm [33].
f  t  1   f  t  P  1    y

3. PROPOSED ALGORITHM
3.1 Overview of the Algorithm(2)

The Nt-dim. vector f is initialized by the source vector, i.e., f(0)=y.
The number of iteration T determines range of relevance diffusionon the manifold. As we will show in the experiments, small T (e.g.,
T=5) is sufficient to obtain high retrieval accuracy. α=[0, 1) is aregularization parameter for relevance diffusion. Since the matrix
P is sparse, computing Equation (2) is quite efficient.

For accurate 3DMR, we propose a novel feature aggregation
algorithm called Diffusion on Manifold (DM) and a novel local 3D
geometric feature Position and Orientation Distribution (POD).
The DM algorithm aggregates a set of local features extracted froma 3D model into a feature vector by using relevance diffusion on amanifold graph of local features (Figure 1d). As a pre-processing
for DM aggregation, the manifold graph is generated from a set oflarge number of training local features. However, the manifold
graph may suffer from burstiness of local features [11], which is aphenomenon where certain types of local features appear more
frequently due to repetitive structures in 3D models or small-scalefeatures that capture too primitive shapes such as flat surfaces.
These “bursty” local features are thought to form dense clusters on
the manifold. In such a case, relevance from bursty local features
would diffuse only within the dense clusters. As a result, theaggregated features would become inaccurate since they are
dominated by relevance values diffused within the bursty clusters.

After relevance diffusion, the Nt-dim. vector f(T) is normalized byusing the similar method to the existing feature aggregation
methods such as FV or VLAD. f(T) is power-normalized by taking
square-root of each element of f(T). Then, the power-normalizedf(T) is normalized by its L2-norm to generate the aggregated feature
vector for the 3D model. Similarity between a pair of two DMaggregated features is computed by using Cosine similarity.

3.2.2 Pre-processing for DM Aggregation
3.2.2.1 Constructing a Manifold Graph
In this section, we describe the method for constructing the
manifold graph P used for DM aggregation. We first randomly subsample the set of Nt (e.g., Nt=250K) training local features from theset of local features extracted from all the 3D models in the
database. A sparse affinity matrix W having size Nt×Nt is generated
by connecting the training local features which are close to eachother in the local feature space. The element W(i, j) of W, which

To alleviate this burstiness problem, we perform weighting ofnodes on the manifold graph. We use Lp-norm IDF weighting
algorithm by Zheng et al. [32] to assign smaller weights to bursty
nodes and to assign larger weights to non-bursty nodes. These node3

Takahiko Furuya and Ryutarou Ohbuchi, Diffusion-on-Manifold Aggregation of Local Features for Shape-based 3D Model Retrieval,
accepted as oral paper, Proc. ACM International Conference on Multimedia Information Retrieval (ICMR) 2015, Shanghai, China.
use the algorithm by Osada et al. [21] for converting a polygonalmodel into an oriented point set. The algorithm randomly anduniformly samples points on the surfaces of the 3D polygonalmodel. Each point is associated with the normal vector of thetriangle on which the point is sampled. In this paper, we sample
Np=3,000 oriented points per 3D model. Oriented point set of the
3D model is scaled to fit a sphere having diameter 1.

indicates the similarity between the local feature i and the local
feature j, is computed by using the following equation;
 exp   d  i , j  /   if i  j and j  kNN  i W i, j   
0otherwise

(3)

where kNN(i) is a set of k nearest neighbor local features of a local
feature i. Using smaller k (e.g., k=5) makes the affinity matrix W
sparser, resulting in faster aggregation. d(i, j) is L1 distance, whose
range is normalized in [0,1], between a pair of two local features. σ
is a scaling parameter for distance. We fix σ to 0.5 in this paper.

3.3.2 Extracting Local Features
Given an oriented point set of the 3D model, we densely extract a
set of POD features from the oriented point set (Figure 2). For eachoriented point, we define a Sphere-Of-Interest (SOI) whose radius
is R. For robustness against scale change of the local 3D shapes, we
use multi-scale SOIs. For each SOI, R is selected randomly from a
range [r-v, r+v]. For example, if we use r=0.5 and v=0.2, R is
randomly selected from a range [0.3, 0.7].

The affinity matrix W is then normalized to generate the probability
transition matrix P as P  D1W where D is a diagonal matrix
whose i-th diagonal element is D  i, i    Nt W  i, j  . The elementj 1

P(i, j) in P indicates the transition probability of relevance from thelocal feature i to the local feature j. P represents the manifold graphof training local features and is used for DM aggregation.

To achieve invariance against 3D rotation of local 3D shapes, wenormalize orientation of the SOIs. For each SOI, we perform PCA
on coordinates of oriented points within the ROI to obtain principle
axes in the rotated coordinate system. We then disambiguate thedirections, or signs, of the principle axes by using the method
similar to SHOT [27] or RoPS [8]. That is, the sign s1 for the firstnprincipal axis e1 is computed as s  sign
 p  p   e where

3.2.2.2 Weighting Nodes on a Manifold Graph
To alleviate the burstiness problem of local features, we compute aweight for each of Nt nodes on the manifold graph P. We use Lpnorm IDF weighting algorithm [32]. It assigns smaller weights tobursty nodes and assigns larger weights to non-bursty nodes.



c



p

i 1

i

1



np is the number of oriented points within the SOI, pi is the positionof i-th oriented point, and p is the center of the SOI. The sign s3 forthe third principal axis e3 is computed in the same manner as s1. The
second principal axis e2 is computed as e2  s3e3  s1e1 where ×

We perform clustering on the set of Nt training local features. For
efficiency, we use the Extremely Randomized Clustering Tree
(ERC-Tree) algorithm [7]. The ERC-Tree algorithm recursively
splits the local feature space until the number of training local
features within the cluster becomes less than Smin. We use Smin=40
in this paper. The weight wc for the cluster c is computed by using
the following equation;

wc  log 1  NM mM fmp,c1denotes cross product of two vectors.
After orientation normalization, we extract a POD feature from therotated SOI. We first compute a bounding box of the SOI, which
may be a cuboid. And we normalize the cuboid bounding box bytransforming it into a cube. We then split the normalized boundingbox into N1×N2×N3 cuboid cells by using a regular grid. N1 (N2, N3)
is the number of cells along the first (second, third) principal axis.
In this paper, we set N1>N2>N3 to make the POD feature compact.
That is, we put finer grids along the principal axis having larger
variance of oriented point distribution. We use N1×N2×N3=4×2×1.

(4)

where NM is the number of 3D models in the database, Mc is a set
of 3D models containing local features that belong to the cluster c,and fm,c denotes the frequency of local features of the 3D model mwithin the cluster c. p is a parameter. The weight wc is assigned to
all the training local features within the cluster c.

For each cell c, we describe the distribution of oriented points. The
position of oriented points within the cell is described by using SV
coding [34]. We count the number of oriented points nc within thecell c and pool displacements of the nc oriented points from thecenter of the cell c. The normals of oriented points within the cellis described by their 3×3 covariance matrix C. These statistics areconcatenated to form a 10-dimensional feature vector for the cell c;f     n , 1 n d , v  where dc is the 3-dim. displacement

Since the ERC-Tree is a randomized clustering algorithm, accuracyof node weights due to a single tree would be insufficient. To obtainreliable node weights, we perform ERC-Tree clustering for Ntree
times and take an average of Ntree weights produced by Ntree trees
to generate final weight wn for the node n. We use Ntree=200 trees
in this paper. The set of Nt weights for all the Nt nodes is used forgenerating the source vector y (Equation 1) for DM aggregation.

3.2.3 Reducing Dimensions of Aggregated Features

c



c



c



c

c



vector pooled within the cell c, and vc is the 6-dim. vector whose
elements are upper triangle elements of the covariance matrix C.

Since DM aggregation tends to produce high-dimensional (e.g.,
Nt=250K dim.) vector, comparing DM-aggregated features is costly.
To accelerate feature comparison, we reduce the dimensionality ofDM-aggregated feature vectors by using Kernel PCA (KPCA) withdot kernel as with [5]. To train the KPCA, we use a set of DMaggregated features extracted from 5,000 randomly selected 3Dmodels in the database. The number of reduced dimensions Nd isless than 100 for most cases, which is much smaller than thedimensionality of the DM-aggregated feature vectors. Comparison
among dimension-reduced features is done by Cosine similarity.

normalize rotation

compute a feature 0.51  0.04    0.11 

3.3 POD Feature for 3D Oriented Point Set

SOI

3.3.1 Generating Oriented Point Set

oriented points

The POD feature extraction algorithm expects 3D model defined asoriented point set. For a 3D model defined as polygonal mesh, we
first convert it into oriented point set by sampling its surfaces. We

cellrotated SOI

POD feature

Figure 2. Extracting POD features from an oriented point set.

4

Takahiko Furuya and Ryutarou Ohbuchi, Diffusion-on-Manifold Aggregation of Local Features for Shape-based 3D Model Retrieval,
accepted as oral paper, Proc. ACM International Conference on Multimedia Information Retrieval (ICMR) 2015, Shanghai, China.
β is a balancing parameter between the frequency and the two
vectors dc and vc. We set β=0.001 in this paper.

4.1.3 Local Features
We compare the POD against five local features; SI [13], LSF [20],and RoPS [8] as local 3D geometric feature, and DSIFT [4, 5] andMO1SIFT [5] as local 2D visual feature. Table 1 summarizes
parameters for DM aggregation used for PSB.

The set of 10-dimensional vectors fc computed for all the N1×N2×N3
cells is concatenated to form a POD feature vector for the SOI. If
we use N1×N2×N3=4×2×1, the POD feature has 80 dim., which ismore compact than other existing local 3D geometric features such
as SI (153 dim.), LSF (625 dim.), FPFH (125 dim.), RoPS (135dim.), and SHOT (320 dim.).

Table 1. Parameters used for DM aggregation (PSB).
POD

k
α
p

4. EXPERIMENTS AND RESULTS
4.1 Experimental Setup
4.1.1 Benchmark Databases

5
0.1
1.5

DSIFT MO1SIFT
5
10
0.01
0.5
1.5
2.5

LSF
5
0.1
1.5

SI

RoPS
5
0.1
2.0

5
0.1
1.5

For local 3D geometric features, i.e., POD, SI, LSF, and RoPS, we
sample Np=3,000 oriented points from surfaces of a 3D model. The
parameters for radius of SOIs are set to r=0.5, v=0.4 for PSB, ESB,
and SH14LC and r=0.1, v=0.1 for SH11NR.

To evaluate accuracy and efficiency of the proposed algorithms, we
use four benchmark databases; the Princeton Shape Benchmark(PSB) [25], the Engineering Shape Benchmark (ESB) [10], the
SHREC 2011 Non-Rigid watertight meshes dataset (SH11NR) [14],and the SHREC 2014 Large-scale Comprehensive 3D shape
retrieval dataset (SH14LC) [15]. Figure 3 shows examples of 3Dmodels as well as the number of 3D models and the number ofcategories contained in the datasets. The PSB and SH14LC containgeneric, rigid 3D shapes. The ESB consists of mechanical CAD
models. The SH11NR has articulated, non-rigid 3D shapes.

For local 2D visual features, i.e., DSIFT and MO1SIFT, a 3Dmodel is rendered from 42 viewpoints spaced uniformly in solid
angle. For DSIFT, We densely and randomly extract a set of 300
SIFT [18] features from a rendered image, resulting in about 13,000
SIFT features per 3D model. To reduce bursty SIFT features, we
use scale-weighted sampling of DSIFT [19]. Specifically, for PSB,
ESB, and SH14LC, we set the scale-weighting parameter W [19] to
0.25, which means more SIFT features are sampled at larger scales
on the image. On the other hand, for SH11NR, we use W=4.0
meaning that more SIFT features are extracted from smaller scales.
MO1SIFT rotates the rendered image to 16 different orientations
and describes each rotated image by a global SIFT feature. Since
we use 42 views for rendering, 42×16=672 MO1SIFT features are
extracted per 3D model.

For all the benchmarks, a 3D model in the database is used as aquery and remaining 3D models are used as retrieval targets. We
use Mean Average Precision (MAP) [%] as an index of retrieval
accuracy. For the SH11NR, to make local feature more robust
against articulation, 3D models are deformed to their canonicalforms by using a method by Jain et al. [9] prior to feature extraction.

4.2 Experimental Results
4.2.1 Effectiveness of DM Aggregation(a) PSB (NM =907, NC =92)(b) ESB (NM =867, NC =45)(c) SH11NR (NM =600, NC =30)(d) SH14LC (NM=8,987, NC=171)

Figure 4 compares retrieval accuracies of the feature aggregation
algorithms for the PSB. We use POD as a local feature. In Figure 4,
we plot MAP scores by varying the number of training local
features Nt for DM aggregation or the number of codewords forother aggregation algorithms. The DM aggregation significantly
outperforms other feature aggregation algorithms. Aggregating
local features considering the non-linear manifold structure of local
features is effective for higher retrieval accuracy.
Table 2, 3, 4, and 5 summarizes retrieval accuracies of feature
aggregation algorithms for the PSB, ESB, SH11NR and SH14LC,
respectively. For benchmarks containing rigid shapes, i.e., PSB,
ESB, SH14LC, the DM aggregation performs the best among theseven feature aggregation algorithms for most local features. DM
aggregation is effective for both local 3D geometric features andlocal 2D visual features.

Figure 3. Benchmark databases used in the experiments(NM : number of 3D models, NC : number of categories).

4.1.2 Feature Aggregation
We compare the DM aggregation algorithm against five existing
feature aggregation methods; BF [2], LL [29], FV [22], VLAD
(VL) [12], and SV [34]. For SV, we use two variants of SV; one iskSV which learns codewords by using k-means clustering, andanother is gSV [5] which learns codewords by GMM clustering andassigns each local feature to multiple codewords according toposterior probabilities that the local features belong to thecodewords.

On the other hand, for SH11NR having non-rigid shapes, the DM
aggregation and its competitors produce similar MAP scores formost cases. We speculate that a key factor behind high accuracy ofDM aggregation is “smoothness” of the manifold graph of local
features. That is, if local features extracted from 3D models
distribute smoothly or continuously in the local feature space,
relevance would diffuses appropriately on the manifold graph oflocal features. For PSB, ESB, and SH14LC, SOIs are highly
overlapped each other since we used SOIs having large radius. In
such a case, local features, whose sampled positions on 3D modelsurfaces are near each other, can also located near in the local
feature space with high possibility. Therefore, the manifold graphof the large-scale local features becomes smooth and DM
aggregation works well for PSB, ESB, and SH14LC. On the other

For BF, we use ERC-Tree clustering to generate about 30K
codewords. For LL, we use k-means to learn 8K codewords. The
number of codewords for FV, VL, kSV, and gSV are determined sothat dimensionality of aggregated feature vector becomes about
300K. For DM, we use Nt=250K training local features. The
aggregated vectors are power-normalized and are then L2normalized. We use Cosine similarity to generate ranking results.

5

Takahiko Furuya and Ryutarou Ohbuchi, Diffusion-on-Manifold Aggregation of Local Features for Shape-based 3D Model Retrieval,
accepted as oral paper, Proc. ACM International Conference on Multimedia Information Retrieval (ICMR) 2015, Shanghai, China.
Table 6 shows effectiveness of node weighting on manifold graphfor the PSB. In Table 6, “without weighting” indicates accuracies
of DM aggregation without node weighting, that is, node weights
wn in Equation 1 is fixed to 1, and “with weighting” shows
accuracies of DM aggregation with node weighting. We canobserve that node weighting improves retrieval accuracy.
Decreasing weights for densely distributed nodes on manifold
graph alleviates the problem of local feature burstiness to generate
more accurate aggregated feature vectors.

hand, the manifold graph for SH11NR might not be smooth since
SOIs with small radius don’t overlap each other. Consequently,
improvement of accuracy due to DM aggregation would be limited
for the SH11NR.
70
65
60

MAP [%]

55

Table 6. Effectiveness of node weighting (MAP [%] for PSB).

50
45

BF-POD
FV-PODkSV-POD
DM-POD

40
35
30
0

LL-POD
VL-POD
gSV-POD

without weightingwith weighting

Table 2. MAP [%] of feature aggregation (PSB).
DSIFT MO1SIFT
54.0
51.9
57.6
56.5
61.7
54.2
60.9
50.6
61.3
49.6
63.8
53.2
64.7
61.1

LSF
33.0
33.8
36.7
38.3
40.1
40.2
41.4

SI
38.1
41.0
44.9
45.5
47.8
49.6
47.3

RoPS
40.5
46.4
46.3
42.7
46.8
48.4
50.4

MAP [%]

POD
47.8
53.1
52.9
52.9
54.8
55.3
60.1

Table 3. MAP [%] of feature aggregation (ESB).
BF
LL
FV
VL
kSV
gSV
DM

POD
52.1
53.0
54.4
49.8
52.4
54.0
57.2

DSIFT MO1SIFT
53.6
49.1
56.0
53.8
57.8
52.2
57.5
47.3
58.7
47.6
59.2
47.9
59.0
57.4

LSF
47.1
50.7
52.8
49.7
51.8
52.6
54.3

SI
48.3
50.9
51.7
46.5
47.1
50.9
54.1

RoPS
47.2
49.1
48.3
44.0
49.9
49.8
51.1

POD
87.3
94.7
96.7
95.9
96.1
96.5
95.8

DSIFT MO1SIFT
94.3
75.7
97.0
82.4
95.6
74.7
96.8
73.1
97.2
74.5
97.0
75.5
96.8
89.2

LSF
87.5
93.8
93.0
95.9
96.0
96.2
95.7

SI
85.8
94.2
95.2
95.8
95.9
95.4
96.9

DSIFT MO1SIFT
38.5
36.5
41.7
39.4
44.1
36.6
43.4
34.9
44.1
35.5
45.5
36.3
44.9
45.5

LSF
30.9
31.4
35.7
34.5
36.2
37.1
38.5

SI
32.7
34.6
38.6
37.1
38.9
40.2
39.3

5
10
15
20
25
Number of iterations T for diffusion

30

Figure 6 plots retrieval accuracies against dimensionality of KPCAprocessed feature for the PSB. For all the three features, i.e., DMPOD, DM-DSIFT, and DM-MO1SIFT, MAP score has a peak ataround 100 dimensions. As the dimensionality of DM-aggregatedfeature is 250K, dimension reduction by KPCA down to about 100
dimensions significantly improve computational efficiency offeature comparison among 3D models.

RoPS
89.0
94.2
94.9
95.2
95.6
93.3
93.8

MAP [%]

POD
39.6
44.7
44.8
44.0
45.2
45.3
50.7

DM-POD
DM-DSIFT
DM-MO1SIFT

Figure 5. Number of iterations for relevance diffusion andretrieval accuracy (PSB).

Table 5. MAP [%] of feature aggregation (SH14LC).
BF
LL
FV
VL
kSV
gSV
DM

70
65
60
55
50
45
40
35
30
0

Table 4. MAP [%] of feature aggregation (SH11NR).
BF
LL
FV
VL
kSV
gSV
DM

DM-DSIFT DM-MO1SIFT
61.4
60.8
64.7
61.1

Figure 5 plots retrieval accuracies against the number of iterations
T for relevance diffusion on the manifold graph. Interestingly, for
DM-POD and DM-DSIFT, accuracies saturate at T=1. While DMMO1SIFT has a slight peak at around T=5. For larger T, retrieval
accuracies are almost unchanged. This is because, presumably, we
use small number of neighbors k for constructing manifold graphs(k=5 for POD and DSIFT, k=10 for MO1SIFT as shown in Table
1) to make feature aggregation efficient. In such a case, themanifold graph would consist of multiple sub-manifolds whichhave almost no connections among them. Relevance diffuses only
within a sub-manifold and thus retrieval accuracy remains
unchanged if we use larger T.

100,000
200,000
300,000
Aggregated feature dimensons

Figure 4. Comparison of feature aggregation methods (PSB).

BF
LL
FV
VL
kSV
gSV
DM

DM-POD
57.8
60.1

RoPS
35.8
39.2
38.6
35.5
38.2
40.1
42.4

DM-POD + KPCA
DM-DSIFT + KPCA
DM-MO1SIFT + KPCA

70
65
60
55
50
45
40
35
30
0

100

DM-POD
DM-DSIFT
DM-MO1SIFT

200 300 400 500 600 700
Reduced feature dimensions Nd

800

900

Figure 6. Reduced dimensions and retrieval accuracy (PSB).

6

Takahiko Furuya and Ryutarou Ohbuchi, Diffusion-on-Manifold Aggregation of Local Features for Shape-based 3D Model Retrieval,
accepted as oral paper, Proc. ACM International Conference on Multimedia Information Retrieval (ICMR) 2015, Shanghai, China.
aggregating the set of local features of the 3D model is still
acceptable for efficient 3DMR.

4.2.2 Effectiveness of POD Feature

MAP [%]

Figure 7 compares retrieval accuracy of POD feature against otherexisting local 3D geometric features, i.e., SI, LSF, and RoPS, for
the PSB. We use DM algorithm to aggregate a set of local features.
In Figure 7, the POD performs the best among the four local
features. Describing distribution of oriented points within SOI by
using SV coding works well for 3DMR.
70
65
60
55
50
45
40
35
30

DM-POD
DM-SI

0

Table 7. Computation time [s] of feature aggregation.
algorithm
VL
gSV
DM

DM-LSF
DM-RoPS

Table 8 shows computation time per query for the SH14LC. In table
8, column “Feat.” is a time for extracting a set of Np=3,000 POD
features from a query 3D model, “Agg.” is a time for aggregating
the set of POD features of the query by using DM, “KPCA” is atime for reducing dimensionality of DM-aggregated feature of thequery from Nt=250K dimensions down to Nd=50 dimensions, and“Sim.” is a time for computing similarities among the query and
8,987 3D models in the database. The proposed DM-POD takes
only about 1.5 [s] to query the SH14LC, which is one of the largest
dataset for 3DMR.

100,000
200,000
300,000
Aggregated feature dimensons

Table 8. Computation time [s] per query for the SH14LC.

Figure 7. Retrieval accuracies for 3D local features (PSB).

algorithm
DM-POD

MAP [%]

Next, we evaluate impact of radius for SOIs on retrieval accuracy.
Figure 8 plots MAP scores, for the PSB, against the parameter r of
SOIs. For single scale POD, we fix another parameter v for SOI
radius to 0. For multi-scale POD, we set v to values described in
Section 4.1.3. We use the DM for aggregating the set of POD
features. For the PSB, accuracies have peaks at around r=0.3 to 0.5
and multi-scale POD shows higher MAP than single-scale POD.
For SH11NR, on the other hand, the peaks exist at r=0.1 and singlescale POD and multi-scale POD performs comparably well. The
articulated models in SH11NR prefer local features extracted fromsmaller radius for robustness against articulation. On the other hand,
the rigid models in the PSB prefers SOIs having more global andmore diverse radius.
100
90
80
70
60
50
40
30
20
0.2

0.4
0.6
Radius r for SOI

0.8

Feat.
0.75

Agg.
0.66

KPCA
0.03

Sim.
0.01

total
1.45

4.2.4 Comparison with Other Retrieval Algorithms
Table 9 compares retrieval accuracy of the proposed DM-POD
algorithm with five state-of-the-art 3DMR algorithms; SV-DSIFT[5], LL-MO1SIFT [5], ZFDR [15] and DBSVC [15]. Of the fivealgorithms, ZFDR combines both 2D visual feature and 3D
geometric feature, and others extract 2D visual feature. We also
evaluate effectiveness of distance metric learning using Manifold
Ranking (MR) [33]. In table 9, “without MR” indicates accuracy
without MR, that is, feature vectors for 3D models are directlycompared using fixed similarity metric such as Cosine similarity togenerate ranking results. And “with MR” shows accuracy using
MR algorithm, that is, ranking results are generated by relevancediffusion on the manifold graph of 3D model features.
The proposed DM-POD shows the highest accuracy when MR isn’t
used. The new local 3D geometric feature combined with the DM
aggregation significantly outperforms other 2D view-based
retrieval algorithms. On the other hand, when MR is applied, DMPOD shows MAP comparable to that of the DBSVC, which is thebest performing algorithm among the SH14LC track entries [15].

single-scale POD (PSB)
multi-scale POD (PSB)
single-scale POD (SH11NR)
multi-scale POD (SH11NR)
0.0

pre-processing
feature aggregation / model
213.4
0.09
2876.2
0.50
167.2
0.66

Table 9. MAP [%] for the SH14LC.

1.0

algorithm
SV-DSIFT [5]
LL-MO1SIFT [5]
ZFDR [15]DBSVC [15]DM-POD (proposed)

Figure 8. Feature scale and retrieval accuracy (PSB).

4.2.3 Efficiency
In this section, we evaluate efficiency of DM aggregation and POD
feature. We use a PC having two Intel Xeon E5-2650V2 CPUs and
256GB DRAM. All the programs run on a single thread.
Table 7 compares computation times for DM, gSV, and VL withrespect to their pre-processing time and feature aggregation time.
We sample Np=3,000 POD features per 3D model. For VL and gSV,
pre-processing means learning 3,000 codewords. For DM, preprocessing includes manifold graph construction and node
weighting on the manifold graph. The DM aggregation shows theshortest pre-processing time as it doesn’t require costly clusteringsuch as k-means or GMM clustering of local features. On the other
hand, the DM is slightly slower than the gSV for aggregating
Np=3,000 POD features per 3D model. However, 0.66 [s] forwithout MRwith MR
46.4
39.9
38.7
44.6
49.5

53.1
46.5
54.1
54.3

5. CONCLUSION AND FUTURE WORKIn this paper, we proposed Diffusion-on-Manifold (DM), a novel
aggregation algorithm for local features, and evaluated it in a shapebased 3D model retrieval (3DMR) setting. The DM exploits, via
manifold graph of features and diffusion distance, non-linear
manifold structure of local features for more accurate aggregation.
Existing feature aggregation algorithms such as Bag-of-Features[2], VLAD [12], Fisher Vector coding [22], Super Vector (SV)coding [34], on the other hand, often fail to capture non-linearity infeature space due to their use of small number of codewords.

7

Takahiko Furuya and Ryutarou Ohbuchi, Diffusion-on-Manifold Aggregation of Local Features for Shape-based 3D Model Retrieval,
accepted as oral paper, Proc. ACM International Conference on Multimedia Information Retrieval (ICMR) 2015, Shanghai, China.
Kosaka, R., Koyanagi, H., Ohbuchi, R., Tatsuma, A. 2014.
Large Scale Comprehensive 3D Shape Retrieval, Proc. EG
3DOR 2014, 131–140.
[16] Liu, L., Wang, L., Liu, X. 2011. In defense of softassignment coding, Proc. ICCV 2011, 2486–2493.
[17] Liu, Y., Zha, H., Qin, H. 2006. Shape Topics: A Compact
Representation and New Algorithms for 3D Partial Shape
Retrieval, Proc. CVPR 2006, 2025–2032.
[18] Lowe, D. G. 2004. Distinctive Image Features from ScaleInvariant Keypoints, IJCV, 60(2), 91–110.
[19] Ohbuchi, R., Furuya, T. 2009. Scale-Weighted Dense Bag ofVisual Features for 3D Model Retrieval from a Partial View
3D Model, Proc. ICCV 2009 workshop on Search in 3D andVideo (S3DV), 63–70.
[20] Ohkita, Y., Ohishi, Y., Furuya, T., Ohbuchi, R. 2012. Nonrigid 3D Model Retrieval Using Set of Local Statistical
Features, Proc. ICME 2012 Workshop on Hot Topics in 3D
Multimedia, 593–598.
[21] Osada, R., Funkhouser, T., Chazelle, B., Dobkin, D. 2002.
Shape Distributions. ACM Trans. on Graphics, 21(4), 807–832.
[22] Perronnin, F., Sánchez, J., Mensink, T. 2010. Improving thefisher kernel for large-scale image classification, Proc.
ECCV 2010, Part IV, 143–156.
[23] Picard, D., Gosselin, P.-H. 2011, Improving Image Similarity
with Vectors of Locally Aggregated Tensors, Proc. ICIP
2011, 669–672.
[24] Rusu, R. B., Blodow, N., Beetz, M. 2009. Fast Point Feature
Histograms (FPFH) for 3D registration, Proc. ICRA 2009,
3212–3217.
[25] Shilane, P., Min, P., Kazhdan, M., Funkhouser, T. 2004. The
Princeton Shape Benchmark, Proc. SMI 2004, 167–178.
[26] Tao, P., Cao, J., Li, S., Liu, X., Liu, L. 2015. Mesh saliency
via ranking unsalient patches in a descriptor space,
Computers & Graphics, Volume 46, 264–274.
[27] Tombari, F., Salti, S., Stefano, L. D. 2010. Unique signaturesof histograms for local surface description, Proc. ECCV2010, 356–369.
[28] Torki, M., Elgammal, A. 2010. Putting Local Features on a
Manifold, Proc. CVPR 2010, 1743-1750.
[29] Wang, J., Yang, J., Yu, K., Lv, F., Huang, T., Gong. Y. 2010,
Locality-constrained Linear Coding for Image Classification,
Proc. CVPR 2010, 3360–3367.
[30] Yang, J., Yu, K., Gong, Y., Huang, T. 2009, Linear spatial
pyramid matching using sparse coding for image
classification, Proc. CVPR 2009, 794–1801.
[31] Yu, K., Zhang, T., Gong, Y. 2009. Nonlinear Learning using
Local Coordinate Coding, Proc. NIPS 2009.
[32] Zheng, L., Wang, S., Liu, Z., Tian, Q. 2013. Lp-norm IDF
for Large Scale Image Search, CVPR 2013, 1626–1633.
[33] Zhou, D., Weston, J., Gretton, A., Bousquet, O., Schölkopf,
B. 2004. Ranking on Data Manifolds, Proc. NIPS 2004.
[34] Zhou, X., Yu, K., Zhang, T., Huang, T.S. 2010. ImageClassification using Super-Vector Coding of Local Image
Descriptors, Proc. ECCV 2010, 141–154.
[35] Zhu, G., Wang, Q., Yuan, Y., Yan, P. 2013. SIFT on
manifold: An intrinsic description, Neurocomputing,113(30), 227–233.

We also proposed Position and Orientation Distribution (POD), a
local 3D geometrical feature for 3D oriented point set. It describes
distribution of oriented points by using SV coding to generatecompact yet accurate local feature.
In experiments, the DM aggregation showed higher accuracy thanthe existing aggregation. The POD feature also outperformed
existing local 3D geometric features for oriented point. As a future
work, we will evaluate effectiveness of the DM aggregation on 2D
image retrieval setting.

6. ACKNOWLEDGMENTS
This research is supported by JSPS Grant-in-Aid for Scientific
Research on Innovative Areas #26120517 and JSPS Grants-in-Aid
for Scientific Research (C) #26330133.

7. REFERENCES[1] Chatfield, K., Lempitsky, V., Vedaldi, A., Zisserman, A.
2011. The devil is in the details: an evaluation of recent
feature encoding methods, British Machine Vision
Conference (BMVC) 2011.
[2] Csurka, G., Dance, C.R., Fan, L., Willamowski, J., Bray, C.
2004. Visual Categorization with Bags of Keypoints, Proc.
ECCV 2004 workshop on Statistical Learning in Computer
Vision, 59–74.
[3] Donoser, M., Bischof, H. 2013. Diffusion Processes for
Retrieval Revisited, Proc. CVPR 2013, 1320–1327.
[4] Furuya T., Ohbuchi. R. 2009, Dense sampling and fastencoding for 3D model retrieval using bag-of-visual features,
Proc. ACM CIVR 2009, Article No. 26.
[5] Furuya, T., Ohbuchi, R. 2014. Fusing Multiple Features for
Shape-based 3D Model Retrieval, Proc. British Machine
Vision Conference (BMVC) 2014.
[6] Gao, Y., Dai, Q. 2014. View-Based 3D Object Retrieval:
Challenges and Approaches, IEEE MultiMedia, 21(3), 52–
57.
[7] Geurts, P., Ernst, D., Wehenkel, L. 2006. Extremely
randomized trees, Machine Learning, 63(1), 3–42.
[8] Guo, Y., Sohel, F., Bennamoun, M., Lu, M., Wan, J. 2013.
Rotational Projection Statistics for 3D Local Surface
Description and Object Recognition, IJCV, 105(1), 63–86.
[9] Jain V., Zhang, H. 2006. Robust 3D Shape Correspondence
in the Spectral Domain, Proc. SMI 2006.
[10] Jayanti, S., Kalyanaraman, Y., Iyer, N., Ramani, K. 2006.
Developing an engineering shape benchmark for CAD
models, Proc CAD, 38(9), 939–953.
[11] Jégou, H., Douze, M., Schmid, C. 2009. On the burstiness ofvisual elements, Proc. CVPR 2009, 1169–1176.
[12] Jégou, H., Douze, M., Schmid, C., P. Perez. 2010.
Aggregating local descriptors into a compact image
representation, Proc. CVPR 2010, 3304–3311.
[13] Johnson, A.E., Hebert, M. 1999. Using spin images for
efficient object recognition in cluttered 3D scenes, PatternAnalysis and Machine Intelligence, 21(5), 433–449.
[14] Lian, Z., Godil, A., Bustos, B., Daoudi, M., Hermans, J.,
Kawamura, S., Kurita, Y., Lavoué, G., Nguyen, H.V.,
Ohbuchi, R., Ohkita, Y., Ohishi, Y., Porikli, F., Reuter, M.,
Sipiran, I., Smeets, D., Suetens, P., Tabia, H., Vandermeulen,
D. 2011. SHREC'11 Track: Shape Retrieval on Non-rigid 3D
Watertight Meshes, Proc. EG 3DOR 2011, 79–88.
[15] Li, B., Lu, Y., Li, C., Godil, A., Schreck, T., Aono, M.,
Chen, Q., Chowdhury, N. K., Fang, B., Furuya, T., Johan, H.,

8