Unsupervised Learning Framework of Interest Point Via Properties Optimization
Pei Yan Yihua Tan Yuan Xiao Yuan Tai Cai Wen
cs.CV

	plain
	
	
	
	
	
	
	Unsupervised Learning Framework of Interest Point Via Properties Optimization
    Pei Yan, Yihua Tan, Yuan Xiao, Yuan Tai, Cai Wen

		National Key Laboratory of Science and Technology on Multi-spectral Information Processing, 
School of Automation, 
		Huazhong University of Science and Technology, Wuhan 430074, China

		yanpei@hust.edu.cn yhtan@hust.edu.cn
		
		
		
		
		
		
		
		
		
		
	
    December 21, 2021
========================================================================================================================================================================================================================================================================================================================


	
	
	
	
		This paper presents an entirely unsupervised interest point training framework by jointly learning detector and descriptor, which takes an image as input and outputs a probability and a description for every image point. The objective of the training framework is formulated as joint probability distribution of the properties of the extracted points. The essential properties are selected as sparsity, repeatability and discriminability which are formulated by the probabilities. To maximize the objective efficiently, latent variable is introduced to represent the probability of that a point satisfies the required properties. Therefore, original maximization can be optimized with Expectation Maximization algorithm (EM). Considering high computation cost of EM on large scale image set, we implement the optimization process with an efficient strategy as Mini-Batch approximation of EM (MBEM). In the experiments both detector and descriptor are instantiated with fully convolutional network which is named as Property Network (PN). The experiments demonstrate that PN outperforms state-of-the-art methods on a number of image matching benchmarks without need of retraining. PN also reveals that the proposed training framework has high flexibility to adapt to diverse types of scenes.
	
	
	
	

§ INTRODUCTION

	
	
	Interest point is a sparse set of image point containing representative information <cit.>. With advantage of computation efficiency and representation stability, interest point has a wide utilization in matching tasks such as simultaneous localization and mapping, image registration and stereo reconstruction. Both computation efficiency and representation stability rely on the location and description of interest point which are produced by detector and descriptor respectively.
	
	All interest point related methods (including detectors and descriptors) can be classified as hand-crafted and learning based methods. Hand-crafted methods define some explicit criteria to extract the points <cit.> and calculate the descriptions while learning based methods are more flexible by training detector and descriptor according to some objectives <cit.>. Both hand-crafted methods and learning based methods make efforts to lead some important properties to interest point which benefit to diverse applications. 
	
	Currently, unsupervised learning is especially feasible to interest point by avoiding the tedious labeling of ground truth points <cit.>. Following the trend, this paper presents an entirely unsupervised training framework by jointly learning detector and descriptor, which takes an image as input and outputs a probability and a description for every image point. Here the probability reflects how likely the point become interest point, and the description is a feature vector used to represent the unique discriminability of an interest point. In essence, the learning model is to approach to the several optimal properties of interest point.  The joint probability distribution of the extracted point properties is used to formulate the objective of our training framework which is maximized to learn detector and descriptor through achieving all desired properties.
	
	The rest of this paper is organized as follows. In section <ref> we discuss the related work. In section <ref>, we formulate our unsupervised learning framework of interest point. Section <ref> introduces latent variable to convert the objective function and explains the Expectation Maximization (EM) optimization algorithm. Section <ref> instantiate detector and descriptor with fully convolutional network, and show some experiments. Finally, in Section <ref> we summarize this paper and list some possible future work.
	
	

§ RELATED WORK

	
	The key idea of hand-crafted interest detector is to define explicit criteria to detect and describe interest point. Harris <cit.>, SIFT <cit.>, SURF <cit.> and KAZE <cit.> extract the points whose gray value have abrupt changes in two-dimensional space or scale space. FAST <cit.> detects the points whose local gray distribution conforms to typical corner patterns. Gradient or gradient-like histograms are widely used in float-value description such as SIFT <cit.>, SURF <cit.>, and Harris <cit.>. Binary descriptors such as BREAF <cit.>,  BRISK <cit.>, FREAK <cit.> select some pixel pairs in neighborhood of interest point to calculate binary feature. Furthermore, scale and orientation estimations <cit.> are normally integrated into above descriptors to improve matching performance under scale and rotation transformations.
	
	TaSK <cit.> and TILDE <cit.> focus on the points extracted by existing detectors (i.e., Forstner <cit.> and SIFT <cit.>) and train detectors to improve repeatability of these points under different illuminations. Quad-network <cit.> approximates repeatability objective with ranking objective and doesn’t rely on any existing detectors. LIFT <cit.> first selects image patches around SIFT point <cit.>, then train descriptor, orientation estimator and fine-tune detector. With known key point locations, SuperPoint <cit.> uses synthetic geometry shapes dataset to pre-trains a weak detector which is further trained in active learning way. Descriptor of SuperPoint is learned by another discriminability objective. In LF-Net <cit.> repeatability objective is approximated by calculating loss of detector between a pair of images, and its discriminability objective impact both descriptions and locations of extracted points, which is similar to LIFT <cit.>. SIPS <cit.> fixes descriptor as existing model such as SIFT <cit.>, and train detector to fit the descriptor. IMIP <cit.> jointly learn detector and descriptor to extract the fixed number of points whose descriptions are limited as one-hot feature vectors.
	
	There are three essential differences between our training framework and existing methods. First, whereas existing methods improve some properties implicitly, we directly maximize the probability that interest point satisfies desired properties. Second, whether a point satisfies the required properties can be thought as latent binary variable so that our training framework can be optimized by Expectation Maximization algorithm. Third, our framework is very flexible to be generalized to any specific application by instantiating different models.
	
	

§ UNSUPERVISED LEARNING FRAMEWORK OF INTEREST POINT

	
	

 §.§ Problem Formulation

	
	Interest point is a sparse set of image point containing representative and discriminative information. In order to integrate both types of properties, a general framework learns detector and descriptor which extract interest point and its description respectively. First we introduce the general process of extracting interest point from a given scene.
	
	For the sake of describing the problem clearly, we give some notations. A actual world scene to acquire image is denoted as I. All possible viewpoint or illumination conditions for taking image are abstracted as transformation set T={T_j | j∈ Z}. Here T_j represents a specific condition, and Z={1,2,…,J} represent all possible conditions. The image acquired from I under condition j is expressed as T_j (I). Each point in I and its corresponding mapped point in T_j (I) are all denoted as o_i. Suppose the entire scene point set is EP={o_i | i=1,2,…,N}. Here N represents the number of scene points. 
	
	In general detector and descriptor take image T_j (I) as input and output interest point and its description respectively. Detector F is defined as a function outputting a probability f_ij for every point o_i in image T_j (I),
	
    f_ij=F(o_i,T_j(I)|θ_F),

	where θ_F are all the parameters of the model of detector, and f_ij reflects how likely o_i becomes an interest point. In practice probability threshold Pt is introduced to obtain a deterministic interest point set. Interest point set of T_j (I) is
	
    IP_j={o_i|f_ij>Pt,o_i∈ EP}.

	And EP-IP_j is named as background point set.
	
	Descriptor D is defined as a function outputting a description vector d_ij for every o_i in image T_j(I), i.e., 
	
    d_ij=D(o_i,T_j(I)|θ_D),

	where θ_D are all the parameters of the model of descriptor, and d_ij can be used to calculate the similarity between this point and other interest point, which is very important to determine the discriminability of an interest point. We always ensure d_ij_2=1 with length normalization. The description set of image T_j (I) is denoted as
	
    DS_j={d_ij|o_i∈ IP_j)}.

	
	The problem of learning based interest point is equivalent to learning the parameters of detector and descriptor. Whereas there is no unique supervised label for interest point and its description, it’s reliable to jointly train detector and descriptor in an unsupervised way which focuses on the essential properties of interest point. 
	Therefore, we propose an unsupervised training framework to optimize the learning model of interest point.
	
	The overview of the framework is shown as Figure <ref>. With images acquired from the same scene, this framework jointly trains detector and descriptor. The training samples can be transformed from a scene image by simulating different illumination and viewpoint changes. Images are fed into detector and descriptor which outputs a probability and a description for every image point. Then, the joint probability of interest point properties is computed through the above two produced information. Finally, a proposed expectation maximization algorithm optimizes the joint probability to find the best model parameters, which means that the outputs of detector and descriptor have  achieved all desired properties.
	
	
	
	
	
	

 §.§ Unsupervised Properties Optimization

	
	The probability that interest point of T_j (I) satisfies given the vth property can be formulated as P_v (IP_j,DS_j). Here v∈{1,2,…,V} and V is the number of all desired properties. Suppose all properties are independent and properties in different images are also independent. The objective maximizing the probability that interest point satisfies desired properties is
	
    argmax_θ_F,θ_D∏_j∏_vP_v(IP_j,DS_j).

	In this paper we make V=3, and sparsity, repeatability and discriminability are selected as essential properties in this paper. The probabilities that an interest point satisfies above properties are denote as sparsity probability, repeatability probability and discriminability probability, whose formulations are introduced in subsequent subsections.
	
	

 §.§ Sparsity Probability

	
	Sparsity is an essential property controlling the limited number of interest points sparsely scatter over image. According to definition sparsity only rely on interest point rather than its description, so its probability can be represented as s(IP_j).
	
	As a specific implementation, we define s_loc (IP_j ) as the probability that interest point set of image T_j (I) is locally sparse. Local sparsity indicates that there is none in the interest point o_i’s neighborhood except itself. U(o_i ) represents the neighborhood of point o_i whose radius is rad. Here rad is a small integral (e.g., rad=4), and U(o_i) doesn’t contain o_i itself. Define · returns the number of elements in a given set. Then the local sparsity probability of o_i is defined as
	
    s_loc(o_i)=
    		 
    			1,     { o_i'|o_i'∈ U(o_i),o_i'∈ IP_j }=0 
    
    			0,     otherwise
    		.

	Suppose the local sparsity probabilities of different interest points is independent, the local sparsity probability of entire interest point IP_j is
	
    s_loc(IP_j)=∏_i∈ IP_js_loc(o_i).

	Except making interest point to be sparse locally, sparsity is also a global properties controlling the number of interest points.  Depending on the number of interest points IP_j, s_num ( IP_j ) represents the probability that number of interest point is reasonable.  We define it as
	
    s_num( IP_j )=
    		 
    			1,      N_min< IP_j <N_max
    
    			0,     otherwise
    		.

	Combining Equation <ref> and <ref>,  the sparsity probability is expressed as
	
    s(IP_j)=s_loc(IP_j)· s_num( IP_j ).

	

 §.§ Repeatability Probability

	
	Repeatability indicates how likely a point can be extracted repeatedly under different viewpoints and illuminations. With this definition repeatability is determined on multiple images acquired from the same scene rather than a single image, so repeatability is defined on I to simplify notations. Denote repeatability probability of point o_i in I as r_i, which represents the probability that o_i is extracted in I, i.e.
	
    r_i=lim_J →∞( 1/J∑_j f_ij).

	where f_ij is the probability output by detector for o_i in image T_j (I) which is defined in Equation <ref>. Then 1-r_i is the probability point o_i belongs to background point set. In the reminder of the paper we ignore the limit of J to simplify notations. Suppose the repeatability of each point is independent on each other, the repeatability probability of IP_j satisfies
	
    r(IP_j)= ∏_o_i∈ IP_jr_i∏_o_i∈ EP-IP_j(1-r_i).

	

 §.§ Discriminability Probability

	
	Discriminability denotes how likely an interest point in one image is more similar to the same point than the other interest points in another image. The similarity between point o_i and o_i' is normally defined as the inner product of their description vector d_ij and d_i'j', which is formulated as
	
    sim_iji'j'=⟨ d_ij, d_i'j'⟩.

	sim_iji'j'∈ [-1,1] because d_ij _2=1. 
	If i=i', sim_iji'j' denotes the similarity between the same point in two images, which is termed as positive pair. Otherwise  when i ≠i', sim_iji'j'  represents the similarity between different points in two images, which is termed as negative pair. 
	
	Denote indicator function as 𝕀 (·) which return 1 if and only if logical operation is true. Define function max returning the maximum of a set and returning –Inf for empty set. Then the discriminability probability of an interest point is
	
    c_i=
    			   1/J(J-1)∑_j ∑_j' ≠ j𝕀(sim_ijij'> . 
       . max({sim_iji'j'|o_i'∈ IP_j',i'≠i })) ) 
    ,

	where o_i' must be an interest point rather than an arbitrary point because only interest point will be considered in matching process. In practice Equation <ref> is too sharp to represent the gap between current and optimal discriminability. So we approximate it with
	
    ĉ_̂î=exp(α (h_i-H)),

	where h_i is the difference between similarity of positive pair and negative pairs, and H is the maximum of h_i. α is a factor controlling the sensitivity of discriminability probability ĉ_i with respect to h_i. The formulation of h_i is
	
    h_i=
    			   1/J(J-1)∑_j ∑_j' ≠ j( min(m_p,sim_ijij')- 
       λ/ IP_j'-o_i ∑_o_i'∈ IP_j'-o_imax(m_n,sim_iji'j') )
    		.

	The formulation of h_i is inspired by descriptor loss in [7]. Here m_p ∈ [-1,1] and m_n ∈ [-1,1] are named as positive margin and negative margin, which can be seemed as the target similarity of positive pair and negative pair. λ is a weight to balance positive pair and negative pair. Because sim_iji'j'∈ [-1,1] the maximum of h_i is H=m_p-λ m_n.
	
	Suppose the discriminability of each point is independent, the discriminability  probability of IP_j is
	
    ĉ(IP_j,DS_j)=∏_o_i ∈ IP_jĉ_i.

	Note in Equation <ref> we don’t concern discriminability between interest point and background point, or discriminability between background points themselves.
	
	

 §.§ Objective of Properties Optimization

	
	During training, we generally have a scene set {I^k| k=1,2,...,K}. Denote interest point set of image T_j (I^k ) as IP_j^k where the descriptor of each point o_i is  d_ij^k. Thus, description set of image T_j (I^k) is DS_kj=(d_ij^k| o_i ∈ IP_kj). Notations r_i, h_i, ĉ_i are also extended as r_i^k, h_i^k, ĉ_i^k to represent repeatability and discriminability for point o_i in scene I^k. Then the objective of properties optimization is
	
    argmax_θ_F,θ_D∏_k∏_j s(IP_j^k)· r(IP_j^k)·ĉ(IP_j^k,DS_j^k),

	which is named as properties objective. Note the description set DS_j^k is ignored in s and r because sparsity and repeatability probability don’t concern description of interest point.
	
	
	

§ OPTIMIZATION AND IMPLEMENTATION

	
	

 §.§ Problem Conversion with Latent Variable

	
	
	It’s hard to straightforward optimize properties objective with conventional gradient based algorithm. In Equation <ref> interest point set IP_j is determined by the probability threshold Pt where the derivative of the logical operation doesn’t exist. Therefore, we introduce latent variable to solve this problem.
	
	Binary latent variable y_i^k is formulated for every point o_i^k in scene I^k and y_i^k=1 if and only if point o_i in any image T_j (I^k) satisfy all desired properties. y^k is the vector whose ith component is y_i^k, and Y is the matrix whose item in kth row and ith col is y_i^k. Point o_i^k is defined as satisfied point if y_i^k=1, and satisfied point set of I^k is {o_i| y_i^k=1}.
	
	In original objective <ref> θ_F and θ_D are optimized to make interest point sets IP_j^k and their descriptions achieve the desired properties, meaning that optimal solution of IP_j^k must be the point set {o_i^k| y_i^k=1}. So replace IP_j^k as {o_i^k| y_i^k=1} in Equation  <ref> won’t change its optimal solution. We first convert the property probability to more compact formulations with y_i^k.
	The number of satisfied points can be calculated with y^k  =∑_i y_i^k . Redefine local sparsity probability for y_i^k as
	
    s_loc(y_i^k)=
    	 
    	y_i^k,     ( ∑_o_i'∈ U(o_i)y_i'^k )=0 
    
    	0,     otherwise
    	.

	For vector y^k we define the result of s_loc (y^k ) is still a vector, whose ith component is s_loc (y_i^k ). 
	
	Redefine sparsity probability of y^k as
	
    s(y^k)=𝕀(N_min<y^k<N_max) ·𝕀(s_l o c(y^k)=y^k)

	The conversion of repeatability probability is straightforward.
	
    r( y^k ) =∏_i(r_i^k)^y_i^k( 1-r_i^k) ^1-y_i^k

	Also we can redefine discriminability probability as
	
    h_i^k(y^k)=
    		   1/J(J-1)∑_j∑_j'  j( min( m_p,sim_iji'j')-
       λ/ y^k ∑_i'  iy_i'^kmax(m_n,sim_iji'j') )

	
    ĉ(y^k, D S_j^k)=∏_iexp(α y_i^k(h_i^k(y^k)-H))

	Note ĉ(y^k, D S_j^k) for any j is same because all images of I^k share the same y^k, so we replace it with ĉ(y^k) in this case. In fact, all property probabilities for different images of I^k are exactly equal by sharing y^k. Then properties objective <ref> can be converted as latent properties objective
	
    argmax_θ_F, θ_D, y∏_k s(y^k) · r(y^k) ·ĉ(y^k).

	According to Equation <ref> sparsity probability s(y^k) is either 1 or 0. To maximize latent properties objective s(y^k) has to be ensured as 1, so that objective <ref> is equivalent to
	
    {[      argmax_θ_F, θ_D, y∏_k r(y^k) ·ĉ(y^k); s . t . N_min<y^k<N_max, s_l o c(y^k)=y^k ].

	With logarithm transformation it’s equivalent to
	
    {[         argmax_θ_F, θ_D, y L=∑_k∑_i l_i^k; s . t . N_min<y^k<N_max, s_l o c(y^k)=y^k ].

	where log-likelihood item
	
    l_i^k=y_i^k log r_i^k+(1-y_i^k) log(1-r_i^k)+α y_i^k(h_i^k(y^k)-H)

	In objective <ref> L is named as log-likelihood function.
	
	
	
	

 §.§ Expectation Maximization of latent properties objective

	
	Objective <ref> can be optimized with Expectation Maximization algorithm (EM). Suppose the optimization contains T times iterations. Let θ_F^0 and θ_D^0 are initial parameters, and the parameters of detector and descriptor are θ_F^t and θ_D^t after t times iterations. 
	
	E-step:
	
	In E-step of iteration t, we need to obtain the expectation of log-likelihood function L with respect to y, which is denoted as E_y(L). To achieve it we first estimate the probability distribution P(y|θ_F^t-1,θ_D^t-1). 
	Because different y^k are independent we only discuss the distribution of y^k directly. The formulation for probability of y^k is
	
    {[  P(y^k | θ_F^t-1, θ_D^t-1)=r(y^k) ·ĉ(y^k); s . t . N_min<y^k<N_max, s_l o c(y^k)=y^k ]..

	Denote all possible y^k form the set
	
    Y^k={ y^k | N_min< y^k <N_max,s_loc (y^k)=y^k }

	Then the probability distribution of y_i^k is
	
    P(y_i^k=1 | θ_F^t-1, θ_D^t-1)=1/Z∑_y^k ∈ Y^k  y_i^k · r(y^k) ·ĉ(y^k),

	where Z=∑_y^k ∈ Y^k r(y^k) ·ĉ(y^k) is the normalization factor. Denote p_i^k=P(y_i^k=1), then the expectation of y_i^k is p_i^k. So the expectation of L is
	
    E_y(L)=∑_k∑_i E_y_i^k(l_i^k)

	Here 
	
    E_y_i^k(l_i^k)=
    			    p_i^k log r_i^k+(1-p_i^k) log(1-r_i^k)+ 
       α E_y^k(y_i^k · h_i^k(y^k)) -α p_i^k H

	
	M-step:
	
	In M-step θ_F^t and θ_D^t are obtained by maximizing the expectation of L. Normally the number of parameters is very huge, so Gradient Ascent algorithm (GA) is selected to achieve M-step in this paper. The gradient is computed as:
	
    ∂ L/∂θ^t-1_F = ∑_k ∑_i p^k_i-r^k_i/J · r^k_i · (1-r^k_i)∑_j ∂ F(i,T_j(I^k))/∂θ ^t-1_F

	
    ∂ L/∂θ^t-1_D = α∑_k ∑_i  ∂ E_y^k(y_i^k · h_i^k(y^k))/∂θ_D^t-1

	
	Similarity sim_iji'j' is computed with inner product which is derivable. So if the detector F and descriptor D are derivable, their parameters θ_F and θ_D can be optimized with GA. Theoretically, from θ_F^t and θ_D^t to θ_F^t+1 and θ_D^t+1 GA may need multiple times updates to achieve convergence. In practice, both Equation <ref> and <ref> lead to very high computational complexity, so we introduce Mini-Batch Approximation as an efficient implementation. All details of of this approximation are outlined in Supplementary Section 1.
	
	
	
	

§ EXPERIMENTAL RESULT

	
	

 §.§ Experiment Setup

	
	Our framework is flexible to be integrated with different models. In this paper we implement detector F and descriptor D as two Fully Convolutional Network <cit.> with Batch Normalization <cit.>, and parameters of their encoders are shared. We name this implementation as Property Network (PN). The architecture of PN is demonstrated in  Supplementary Section 2.
	
	In this paper we select MS-COCO 2014 <cit.> as training dataset, which comprises of more than 80 thousand images. We treat each single image as a scene which is transformed to the training images with different illuminations and viewpoints through simulated transformations.  The details of simulated transformations are outlined in Supplementary Section  3. 
	
	To adapt to different kinds of illumination and viewpoints changes, we train three different PN models corresponding to three transform simulation conditions. PN-i-64 model is trained with the condition of sharp illumination change and medium viewpoint change. PN-v-64 is trained with the condition of large range viewpoint change and medium illumination change. PN-128 is trained with the condition of sharp illumination change and large range viewpoint change. The architectures of above three networks are identical except the length of their description vector. The length of description vector of PN-i-64, PN-v-64 and PN-128 are 64, 64, 128 respectively.
	
	The configurations of all hyperparameters of PN are as below. During training we resize all images to 320×240 and stack different images to mini-batches (but the size of testing image can be arbitrary). In every iteration two scenes are randomly transformed ten times respectively, so there are in fact twenty images in a mini-batch (i.e., B=2 and J=10). In entire training we fix number range N_min=200, N_max=400, and non-maximum suppression radius rad=4 pixels. For discriminability, negative pair weight λ=10/N_max , positive margin m_p=1, negative margin m_n=0.2 and discriminability weight α=1. We always stop training after two epochs. The FCN model is implemented with and solved by Adam optimizer <cit.> with default parameters (lr=0.001 and β=(0.9,0.999)).
	
	

 §.§ Performance Comparison

	
	
	Three datasets HPatches <cit.>, Webcam <cit.>, Oxford <cit.> are used to evaluate the performance. HPatches is divided into illumination changed subset and viewpoint changed subset, which are denoted as HP-i and HP-v respectively. Webcam contains sharp illumination changes but no viewpoint change. In this paper no method is trained on Webcam, so training set and testing set of Webcam are both available to evaluate performance, which is denoted as W-train and W-test. Oxford comprises of both illumination and viewpoint changes, but we don’t split it because the number of images in Oxford is small.
	
	Two used metrics Matching Score and Homography Estimation are identical to that used in <cit.> (see more details in Supplementary Section 4). Higher Matching Score and Homography Estimation are better. Briefly, Matching Score measures the ratio of the recovered ground truth correspondences over the number of points extracted by the detector in the shared viewpoint region. Homography Estimation measures the ability of an algorithm to estimate the homography that relates a pair of images by comparing it to ground truth homography. Similar to <cit.>, we fix correct distance ϵ=3 pixels, and use RANSAC method implemented by Opencv Toolbox to estimate homography. In this paper all evaluations are performed on image size 640×480. To be fair we use the recommended hyperparameters for every method, except the maximum number of extracted interest points.  We keep no more than 1000 interest points for 640×480 images. 
	
	The performances of different methods are shown in Table <ref> and Table <ref>.
	
	In summary, our Property Networks outperform the other methods on viewpoint changed subset of HPatches and Oxford, and achieve competitive performance on others according to Matching Score. With Homography Estimation, our methods outperform other methods on HPatches and Webcam and achieve competitive performance on Oxford. As the state-of-the-art algorithms, SIFT has low performance on illumination changed datasets and SuperPoint has no superiority on viewpoint changed datasets, but our Property Networks present much more stable results on all above datasets. 
	Furthermore, it’s interesting that SuperPoint has slight advantage on Webcam according to Matching Score but has large gap comparing to our PN-i-64 according to Homography Estimation. This is because Matching Score is a basic metric only reflecting the maximum number of points which are likely helpful to subsequent tasks. However, the space distribution of matched points is as important as the number of matched points for the final tasks such as homography estimation and camera pose estimation, which can’t be revealed by Matching Score. Sparsity is an essential property making the space distribution of interest point suit for these tasks. By improving sparsity and other properties jointly PN-i-64 achieve much higher performance for Homography Estimation on illumination changed dataset.
	
	Figure <ref> demonstrates some visualization details about it. Whereas SuperPoint and SIFT achieve the best performance in some specific scenes, PN-i-64 model give reliable results under both illumination and viewpoint changes. Furthermore, above visualization results of PN-i-64 and SuperPoint explains the difference between Matching Score and Homography Estimation. 
	More results can be found in Supplementary Section 5.
	
	
	
	

 §.§ Performance Analysis of Three PM Models

	
	
	It’s an ultimate goal to learning an interest point model that can be generalized to all application scenes. From the viewpoint of actual application, it’s practical to learning an interest model to cope with specific conditions. The results of our three models confirm this opinion, which is intuitively demonstrated in Figure <ref>.
	
	Though PN-128 has the highest potentials with largest description vector length, in our training it can’t converge well facing both sharp illumination changes and large range viewpoint changes, which makes it have no superiority for either illumination or viewpoint changes. From another perspective, if there are both sharp illumination changes and large range viewpoint changes in given application, PN-128 should be a reasonable tradeoff with limited learning ability. Note PN-128 achieves similar or better results comparing with SuperPoint whose length of description vector is 256.
	
	By focusing on sharp illumination changes and medium viewpoint changes, PN-i-64 achieves much better performance on illumination changed datasets, which demonstrates the flexibility of our training framework. What we need is only feeding the model corresponding images for different application scenes, without adjusting the architecture and objective anymore.
	PN-v-64 outperform PN-128 under viewpoint changes, but the superiority is relatively small. One reason is conventional convolution neural network such as Fully Convolutional Network can only achieve limited rotation invariant with convolution and pooling operation <cit.>. Explicitly estimating the orientation of interest point is a one of reliable ways to improve rotation invariant <cit.>.
	
	
	
	
	

§ CONCLUSION

	
	This paper proposes an entirely unsupervised training framework by maximizing the sparsity, repeatability and discriminability probability of interest point and its description. With Expectation Maximization algorithm and mini-batch approximation this framework can be optimized efficiently. As an implementation based on Fully Convolutional Network, Property Network outperforms state-of-the-art algorithms on a number of image matching dataset, which demonstrates the effectiveness and flexibility of our training framework. Future work will investigate more well-designed architecture of detector and descriptor which can reach more potential of this framework. Furthermore, our framework can be integrated with more properties to improve model performance in diverse applications, and some supervised information can also be formulated as properties which bring semantics to interest point.
	
	ieee
		
	

	

§ SUPPLEMENTARY MATERIAL

	In this supplementary material we give more details of our implementation and experimential results. In section <ref> we implement the Expectation Maximization (EM) process with an efficient strategy called Mini-Batch approximation of EM (MBEM). Section <ref> demonstrates the architecture of our Property Network. Section <ref> outlines details of our simulation for illumination and viewpoint changes. In Section <ref> we introduce performance metrics used in experiments and more experimential results are demonstrated in Section <ref>.
	
	
	
	

§ MINI-BATCH APPROXIMATION OF EXPECTATION MAXMIZATION

	
	

 §.§ Problems of Original Expectation Maxmization

	
	Before discuss the difficultes of optimizing our objective with original Expectation Maxmization algorithm (EM), we first overview the objective of our training framework, which is formulated as
	
    {[ max_θ_F, θ_D, y^1,.,y^K P(θ_F, θ_D,y^1,.,y^K)=∏_k r(y^k) ·ĉ(y^k);                        s . t . N_min<y^k<N_max, s_l o c(y^k)=y^k ].,

	where θ_F and θ_D are parameters of detector and descriptor respectively, and y^1,...,y^K are latent variables. Each y^k is correspond to a scene I^k, and y^k is a vector whose component y_i^k is a binary random variable representing whether point o_i^k satisfies desired properties. Point o_i^k is defined as satisfied point if y_i^k=1. Note y^k need to be optimized because we don't known which points can become satisfied points. To be convenient for optimization we also define log-likelihood function
	
    L=log(P(θ_F, θ_D,y^1,.,y^K)).

	
	Because each y^k is independent in objective (<ref>) , we drop the superscript in the following parts of the supplement for the brief expression. Without ambiguity we always use L to represent log-likelihood function no matter whether the k is ignored or not.
	
	Properties considered in objective (<ref>) are sparsity, repeatability and discriminability. The repeatability and discriminability probability is
	
    r(y) =∏_i(r_i)^y_i( 1-r_i) ^1-y_i,

	
    ĉ(y)=∏_i(ĉ_i(y))^y_i,

	where r_i and ĉ_i(y) are probabilities that o_i satisfies repeatability and discriminability respectively. The formulations of r_i and ĉ_i(y) can be found in main text. In this section we only need know r_i don't rely on y, but ĉ_i(y) depend on vector y.
	
	The constraint in objective (<ref>) is named as sparsity constraint.  y=∑_iy_i represents the number of satisfied points, and [N_min,N_max] are reasonable number range of interest points in a single scene. Local sparsity constraint s_loc is defined as
	
    s_loc(y_i)=
    	 
    	y_i,     ( ∑_o_i'∈ U(o_i)y_i')=0 
    
    	0,     otherwise
    	,

	where U(o_i ) represents the neighborhood of point o_i whose radius is rad. Here rad is a small integral (e.g., rad=4), and U(o_i) doesn’t contain o_i itself.
	
	Theoretically, objective (<ref>) can be optimized with EM. Normally EM contains T iterations, and an Expectation step (E-step) and a Maximization step (M-step) are conducted in each iteration. Let θ_F^0 and θ_D^0 are initial parameters, and the parameters of detector and descriptor are θ_F^t and θ_D^t after t iterations. In each E-step we need obtain p_i=P(y_i=1 | θ_F^t-1, θ_D^t-1) and compute the expectation of L with respect to y which is denoted as E_y(L). In each M-step we need update θ_F^t-1 and θ_D^t-1 to θ_F^t and θ_D^t by maximizing E_y(L).
	
	We first formulate p_i theoretically. Without sparsity constraint the sample space of y can be represented as {y^(m)| m=1,2,...,2^N} where N is the number of all the points in given scene. That's because the length of vector y is N and each component of y is a binary variable. But with sparsity constraint this sample space should be more narrow. Denote the sample space of y as Y which can be formulated as
	
    Y=
    	   {y^(m)| m=1,2,...,2^N,. 
       . N_min<y^(m)<N_max, s_l o c(y^(m))=y^(m)}.

	
	The sparsity constraint in objective (<ref>) can be ignored if we ensure y∈ Y. So when y∈ Y the distribution of y can be reformulate as
	
    P(y | θ_F^t-1, θ_D^t-1)
    	   = r(y) ·ĉ(y),
       =∏_ir_i^y_i·(1-r_i)^1-y_i·ĉ_i(y)^y_i
 
	Note we obtain (<ref>) by substituting (<ref>) and (<ref>) into (<ref>).
	
	Because the distribution p_i is a marginal distribution of P(y | θ_F^t-1, θ_D^t-1), so we obtain
	
    p_i
    	   =P(y_i=1 | θ_F^t-1, θ_D^t-1)
       =1/Z∑_y^(m)∈ Y,y_i=1 P(y^(m) | θ_F^t-1, θ_D^t-1)
       =1/Z∑_y^(m)∈ Y_y_i=1 r(y^(m)) ·ĉ(y^(m)) 
       =1/Z∑_y^(m)∈ Y_y_i=1∏_ir_i^y_i^(m)·(1-r_i)^1-y_i^(m)·ĉ_i(y^(m))^y_i^(m)
 
	where Z=∑_y^(m)∈ Y r(y^(m)) ·ĉ(y^(m)) is the normalization factor, and Y_y_i=1={y^(m)| y^(m)_i=1, y^(m)∈ Y}. By combining <ref>, <ref> and <ref>, the expectation of L with respect to y is
	
    E_y(L)=
    	   ∑_i(p_i log r_i+(1-p_i) log(1-r_i)+
       E_y(y_ilog(ĉ_i(y)))),

	Note we drop all the k of L and only considering single y. Here ĉ_i(y) depend on y so we can't straightforward simplify the expectation E_y(y_ilog(ĉ_i(y))) in (<ref>).
	
	The E_y(L) need be maximized in M-step, which can be achieved with Gradient Ascent algorithm (GA). According to definitions in main text, r_i depends on θ_F, and ĉ_i(y) depends on θ_D. The expressions of partial derivative are
	
    ∂ L/∂θ^t-1_F = ∑_i p_i-r_i/r_i · (1-r_i)·∂ r_i/∂θ ^t-1_F,

	
    ∂ L/∂θ^t-1_D = ∑_i ∂ E_y(y_ilog(ĉ_i(y))/∂θ_D^t-1.

	Here ∂ r_i/∂θ_F in (<ref>) can be computed directly, but the partial derivative in (<ref>) is more complex.
	
	Above E-step and M-step lead to high computational complexity from three aspects. 1) It’s inefficient to make GA achieve convergence on entire image set in each M-step (note (<ref>) and (<ref>) are defined for one scene but in fact we need compute them for all scenes). 2) Directly computing p_i with (<ref>) need traverse y^(m) through entire Y, which leads to very high computational complexity. 3) Even with known p_i, computing (<ref>) still need traverse y^(m) through entire Y.
	
	In the subsequent subsection, we solve above three problems with some approximations of original EM. The entire approximate algorithm is named as Mini-Batch approximation of Expectation Maxmization (MBEM).
	
	

 §.§ Considering Mini-Batch in Each Iteration

	
	Because it's too slow to make GA achieve convergence on entire image set in each M-step, in each iteration we only focus on a small subset of training set, which is normally named as mini-batch.
	
	Here is more details. First we select size of mini-batch as BS (in our training we fix BS=2). Construct B_t={b | t· BS+1≤ b≤ (t+1)BS} in iteration t and select mini-batch {I^k | k ∈ B_t} to conduct E-step and M-step. Furthermore, each M-step only updates parameters once rather than multiple times until converge, because it’s not necessary to achieve convergence on a mini-batch which increases the risk of overfitting.
	
	

 §.§ Efficient Approximation for Distribution of Latent Variable

	
	Though distribution p_i can be computed with (<ref>), but this computation need traverse y^(m) through entire sample space Y which leads to very high computational complexity. There are two obstacles to simplify (<ref>).
	
		
  * The definition of sample space Y is integrated with sparsity constraint. It's easy to check whether a specific y^(m) satisfies sparsity constraint or not, but it's difficult to straightforward obtain all satisfied y^(m).
		
  * To solve distribution p_i, we must consider all points because ĉ_i depends on vector y.
	
	
	Whereas it's hard to solve (<ref>) precisely, this subsection introduce an efficient strategy to approximate p_i which comprises of three step. First, we approximate sample space Y with Ŷ which can be obtained efficiently. Second, the ĉ_i(y) is approximated with ĉ_i(ŷ) where ŷ is a constant for given θ_F. Third, we simplify the approximate formulation of p_i.
	
	

  §.§.§ Efficient Approximation of Sample Space

	
	We first approximate Y with Ŷ to simplify the sparsity constraint. A reasonable Ŷ should have two characteristics. First, all y∈Ŷ satisfy sparsity constraint. That's mean Ŷ can only be a subset of Y, i.e., we can only remove some elements from Y in order to obtain Ŷ. Second, replace Y with Ŷ won't change p_i significantly. That's mean elements in Ŷ should have greater contributions to p_i comparing to elements in Y-Ŷ.
	
	The first characteristic can be achieved by setting y_i'=0 for all o_i'∈ U(o_i) if we want to setting y_i=1. The remaining question is how to select point o_i to set y_i=1. According to (<ref>), the larger r_i and ĉ_i(y) point o_i has, the greater contribution y_i=1 can make. Considering r_i can be obtained directly but ĉ_i(y) depend on y, and we have assume r_i and ĉ_i(y) are independent in main text, so we can select point o_i if r_i is a local maximum and set y_i=1.
	According to above analyses we define
	
    ŷ_i=
    	
    	1,     r_i > max(r_i'|o_i'∈ U(o_i))
    
    	0,     otherwise
    	,

	where max(·) return the maximum of given set. Define ŷ is a vector whose ith component is ŷ_i. Since the local maximum points that are deducted from the current model require the corresponding points in the sample y^m must be also local maximum as ŷ, the sample space Y can be further reduced as
	
    Ŷ={ y^(m)| N_min< y^(m)<N_max,  ∀ i, y^(m)_i≤ŷ_i } .

	Note any y ∈Ŷ satifies s_loc (y)=y becuase s_loc(ŷ)=ŷ, this significantly benefit computation efficiency by avoid checking this constraint.
	
	

  §.§.§ Efficient Approximation of ĉ_i(y)

	
	The second problem is ĉ_i depends on vector y which make us have to consider all points when solve p_i. To reduce the computational complexity, we give an approximation of ĉ_i(y) with ĉ_i(ŷ).
	We first make an assumption called
	
	Discriminability Consistency:
	
	Given two feasible samples y^(1),y^(2)∈Ŷ which satisfy {o_i | y^(1)_i=1 }⊆{o_i | y^(2)_i=1 }, discriminability consistency 
		assumes ∀ i, i' ∈{o_i | y^(1)_i=1 }, if ĉ_i(y^(1))<ĉ_i'(y^(1)), then ĉ_i(y^(2))<ĉ_i'(y^(2)).
	
	Discriminability consistency assumes if the discriminability of point o_i is larger than o_i' in satisfied point set {o_i | y_i=1 }, then inserting some interest points into this satisfied point set may not change their relative order. According to the definition ŷ which is deduced from the local maximum r computed with current model parameters, ŷ must be the super set of all the other genuine interest point sets.
	
	If we view  any sample y as y^(1) and ŷ as y^(2), then the discriminability consistency will be satisfied. Our objective of optimization is to maintain the best descriminability of each interest point so that we can use ŷ_i to replace y_i. Therefore,it's reasonable to select ŷ to replace any sample  y to guarantee discriminability consistency. Furthermore, the feasible samples are also replaced by ŷ such that we only concern about the computation of ŷ. With ŷ we can obtain 
	
    ĉ(ŷ)=∏_iexp(αŷ_i(h_i(ŷ)-H)),

	where h_i is the function computing disciminability for given point o_i, and its formulation can be found in main text.
	
	

  §.§.§ Simplify the Formulation of p_i

	
	First we give the formulation of p_i combining above two approximations. Because ŷ_i is a constant with given r_i, we denote ĉ_i(ŷ) as c̃_i for the simplicity. Then (<ref>) can be approximated with
	
    p_i≈1/Z∑_y^(m)∈Ŷ_y_i=1∏_i'r_i'^y^(m)_i'·(1-r_i')^1-y^(m)_i'·c̃_i'^y^(m)_i'.

	Here Z is still the normalization factor and we don't show its formulation repeatedly, and Y_y_i=1 in (<ref>) is changed to Ŷ_y_i=1={y^(m)| y_i=1, y^(m)∈Ŷ}.
	
	In order to ensure computation efficieny we need avoid the summation over y^(m), the product over i' and the computation of normalization factor Z in (<ref>). In this subsection we intorduce our solution to achieve it.
	
	We make z_i^(m)=r_i^y_i^(m)·(1-r_i)^1-y_i^(m)·c̃_i^y_i^(m) for the simplicity , then (<ref>) can be reformulated as
	
    p_i
    	   =1/Z∑_y^(m)∈Ŷ_y_i=1 r_i·c̃_i∏_i' iz_i'^(m)
       =1/Zr_i·c̃_i∑_y^(m)∈Ŷ_y_i=1∏_i' iz_i'^(m).

	Then we formulate 1-p_i according to (<ref>), i.e.,
	
    1-p_i
    	   =P(y_i=0)
       =1/Z∑_y^(m)∈Ŷ_y_i=0 (1-r_i)∏_i' iz_i'^(m)
       =1/Z(1-r_i)∑_y^(m)∈Ŷ_y_i=0∏_i' iz_i'^(m),

	where Ŷ_y_i=0={y^(m)| y_i=0, y^(m)∈Ŷ}.
	
	Make · return the number of elements in given set, then (<ref>) and (<ref>) can be converted to
	
    p_i=1/Zr_i·c̃_i·Ŷ_y_i=1· avg_1(∏_i' iz_i'),

	
    1-p_i=1/Z(1-r_i)·Ŷ_y_i=0· avg_0(∏_i' iz_i').

	where 
	
    avg_1(∏_i' iz_i')=1/Ŷ_y_i=1∑_y^(m)∈Ŷ_y_i=1∏_i' iz_i'^(m),

	
    avg_0(∏_i' iz_i')=1/Ŷ_y_i=0∑_y^(m)∈Ŷ_y_i=0∏_i' iz_i'^(m).

	
	Whether y_i=1 or y_i=0 only slightly change ∏_i' iz_i', and this difference can be ignored after average them over Ŷ_y_i=1 and Ŷ_y_i=0, so avg(Ŷ_y_i=1) ≈ avg(Ŷ_y_i=0). With (<ref>) and (<ref>) we can obtain
	
    p_i
    	   =p_i/p_i+(1-p_i)
       ≈r_i·c̃_i·Ŷ_y_i=1/r_i·c̃_i·Ŷ_y_i=1+(1-r_i)·Ŷ_y_i=0

	Note r_i and c̃_i can be computed efficiently without interaction to Ŷ. So only Ŷ_y_i=1 and Ŷ_y_i=0 need to be determined in (<ref>).
	
	We first compute Ŷ with (<ref>). Each y in Ŷ corresponds to a set of satisfied points which is select from { o_i|ŷ_i=1}. This selection can be splited into two steps. First, we determine the number of satisfied points, i.e., select n ∈ [N_min,N_max]. Second, we choose n point from { o_i|ŷ_i=1}, and the corresponding y is obtained. According to above two step Ŷ can be calculated with
	
    Ŷ   =C_ŷ^n
       =∑_n=N_min^N_maxŷ!/n!(ŷ-n)!,

	where C_ŷ^n indicate the number of combinations.
	
	In computation of Ŷ_y_i=1, we can only choose y-1 satisfied points because we have select o_i as a satisfied point. So
	
    Ŷ_y_i=1   =C_ŷ-1^n-1
       =∑_n=N_min-1^N_max-1(ŷ-1)!/n!(ŷ-n)!.

	Then 
	
    Ŷ_y_i=0=Ŷ-Ŷ_y_i=1.

	
	Substituting (<ref>), (<ref>) and definitions of r_i and c_i (details can be found in main text) into (<ref>), the approximate p_i can be computed efficiently.
	
	

 §.§ Efficient Approximation for Partial Derivative of Description Parameters

	
	In Subsection <ref> we can obtain distribution p_i efficiently.
	But even with known p_i computing E_y(y_ilog(ĉ_i(y)) and its partial derivative in (<ref>) still need traverse y through entire Y because ĉ_i depends on y. Fortunately, in Subsection <ref> we have solve this problem by approximating y with ŷ under the assumption of discriminability consistency. Because ŷ is independent of y, so
	
    E_y(y_ilog(ĉ_i(ŷ))
    	   =log(ĉ_i(ŷ))E_y(y_i)
       =p_ilog(c̃_i),

	where c̃_i=ĉ_i(ŷ) is defined in Subsection <ref>. Substituting (<ref>) into (<ref>) we obtain
	
    ∂ L/∂θ^t-1_D   = ∑_i ∂ E_y(y_ilog(ĉ_i(y))/∂θ_D^t-1
       ≈∑_i p_i∂log(c̃_i)/∂θ_D^t-1
       = ∑_i p_i/c̃_i∂c̃_i/∂θ_D^t-1.

	According to the definition of discriminability in main text, ∂c̃_i/∂θ_D^t-1 can be computed directly. So the entire Gradient Ascent algorithm can be performed.
	
	

§ ARCHITECTURE OF PROPERTY NETWORK

	
	
	
	Property Network (PN) is a specific implementation of our training framework. In PN, both detector and descriptor are implemented with Fully Convolutional Network <cit.>, whose architecture is inspired by <cit.>. Figure <ref> shows the architecture of PN briefly. Detector comprises of encoder and detection decoder, and descriptor comprises of encoder and description decoder. Both detector and descriptor share the same encoder. Table <ref> denomstrate more details of architecture of PN.
	
	In Table <ref>, We ignore all Batch Normalization and ReLU layers, which follow every convolution layer except the final convolution layer in Conv5 and Conv4-2. “Concatenate” operation concatenates current feature map and given feature map into a new feature map. "trans-conv" indicate transposed convolution <cit.>. p=1 means padding is 1, and s=1 indicates the strider of sliding window is 1. des_len is the length of description vector. As introduced in main text, we make des_len=64 in PN-i-64 and PN-v-64, and des_len=128 in PN-128.
	
	
	
	

§ SIMULATIONS FOR VIEWPOINT AND ILLUMINATION CHANGES

	
	In this section we outline our simulations for illumination and viewpoint changes. For illumination simulation we randomly select the transformations to change pixel value. For viewpoint simulation we randomly generate homography matrices used to perform homography transformation.  Figure <ref> shows several examples of our simulated results.
	We first introduce more details of simulations for illunimation changes. Our slmulations comprise of seven kinds of transformations. In the discussion below we assume the original image have three color channels with 256 gray levels.
	
		
  * Image Blur. Apply Gaussian blur, average blur or median bluron image.
		
  * Channels Shuffle. Permute the order of the color channels of image.
		
  * Contrast Normalization. Change the contrast in images by moving pixel values away or closer to 128.
		
  * Grayscale. Convert images to grayscale and mixe with the original image with a random weight.
		
  * Invert all pixels in given image, i.e. set them to 255-original_pixel_value.
		
  * Salt and Pepper noise. Randomly replace some pixels with very white or black colors.
		
  * Shadow. Randomly insert some dark shapes into image.
	
	
	In main text we have mentioned different PN models are trained with different level of simulations. PN-v-64 is trained with only Image Blur, Contrast Normalization and Shadow, and PN-i-64 and PN-128 are trained with all kind of simulations for illumination changes.
	
	In simulations for viewpoint changes, we randomly generate homography matrices and use them to perform homography transformation. Different PN models are also trained with different level of homography transformation. In training of PN-i-64, we ensure the rotation angle of homography transformation to be less than 45^∘, and for PN-v-64 and PN-128 we don't conduct restriction of rotation angle.
	
	All above simulations are online, i.e., in each training iteration we randomly perform above transformations on current image mini-bacth, and all transformed images will be discarded after this iteration.
	
	

§ PERFORMANCE METRICS

	
	Two metrics used in our experiment are Matching Score and Homography Estimation, which are identical to that used in <cit.>. In this section we introduce their definitions.
	

 §.§ Matching Score

	Matching Score measures the overall performance of  interest point detector and descriptor. It measures the ratio of ground truth correspondences that can be recovered by detector and descriptor over the number of interest point proposed by the pipeline in the shared viewpoint
	region. Suppose the detector extracts IP_1 points from image I_1 and IP_2 points from image I_2 in the shared viewpoint
	region of I_1 and I_2. With descriptions of extracted points and a specific matching strategy, suppose corr_12 points in I_1 are correctly matched to points in I_2, and corr_21 points in I_2 are correctly matched to points in I_1. Then Matching Score for this pair of image is
	
    M-score=1/2(corr_12/ IP_1+corr_21/ IP_2).

	Same to <cit.>, we use two-way nearest neighbor matching as the matching strategy. And Matching Score for entire dataset is obtained by averaging M-score over all image pairs.
	
	

 §.§ Homography Estimation

	Homography Estimation measures the ability of an algorithm to estimate the homography that relates a pair of images by comparing it to ground truth homography. Denote the estimated homography as Ĥ and ground truth homography as H. Be identical to <cit.>, Homography Estimation is defined with the error of four corners of one image onto the other. Denote the four corners of
	the first image as c_1, c_2, c_3 and c_4. We then apply the ground truth H to get the ground truth corners in the second image which is denote as c_1', c_2', c_3' and c_4', and the estimated homography Ĥ to get the estimated corners in the second image which is denoted as ĉ_̂1̂', ĉ_̂2̂', ĉ_̂3̂' and ĉ_̂4̂'. And the error of homography is defined as
	
    Homo_error=1/4∑_i=1^4ĉ_̂î'-c_i'.

	And the measure of Homography Estimation is defined as
	
    HE=𝕀(Homo_error≤ϵ).

	Here 𝕀 is the indicator function, and ϵ is the threshold to judge the correctness of homography. Same to <cit.> we set ϵ=3 in our experiment.
	
	

§ MORE EXPERIMENTIAL RESULTS

	
	
	
	In main text we have demostrated results on image size 640× 480. In this subsection we give results of different methods on image size 320× 240. All experiment configures are same as that for image size 640× 480 except the maximum number of extracted interest points.  We keep no more than 300 interest points for 320×240 images. Table <ref> and <ref> demostrate Matching Score and Homography Estimation of different results. And Figure <ref>  shows some matching results with image size 320× 240. Here our PN-i-64, PN-v-64 and PN-128 are identical to what in main text.
	
	
	
	
	
	
